Lightchaser

三人行,必有我师

三阶力常数矩阵

需要用到thirdorder程序。

使用方法:

  1. 准备好POSCAR,INCAR,KPOINTS,POTCAR和提交任务脚本(以sub.sh为例);
    其中POSCAR与二阶力常数矩阵求解时的POSCAR-unit一致,KPOINTS也与二阶一致。INCAR文件与二阶力常数求解时的有限位移法一致,其他不变。完成后运行thirdorder实现扩胞。

~/安装目录/thirdorder_vasp.py sow 2 2 2 -12
2 2 2是扩胞倍数,要和二阶力常数矩阵保持一致或者略小。-12指的是考虑多少个近邻原子,需要测试,但是如果计算资源紧缺的时候取个-10就行了,但有风险;正值是指截断半径。之后生成成功以后会出现一堆3RD开头的文件,那些都是POSCAR,一般有几百上千个。

  1. 参考phonopy计算声子谱的有限位移法,需要把这些POSCAR单独计算。也就是有多少POSCAR,建立多少文件夹,分别进行多少次计算。
  2. find job* -name vasprun.xml|sort -n|/ thirdorder_vasp.py reap a b c -d a b c和-d与前面一致。该命令可获得FORCE_CONSTANT_3RD也即是三阶力常数矩阵。

脚本方法

  • 测试-d参数(秦光照老师)

    {.line-numbers}
    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
    #!/bin/bash
    #
    # Construct the 3nd files based on the finished calculations
    #
    # Written BY QIN, GuangZhao <qin.phys@gmail.com>
    # 2016-11-22

    DIR="/vol-th/home/xhaiyan/WangN/Au6S2/third241-14"

    if [ ! -d "$DIR" ];then
    echo Please specify the source directory.
    exit
    fi

    for NAME in 3RD.POSCAR.*
    do
    NUM="${NAME#3RD.POSCAR.}" # the number of the 3nd POSCAR file
    CODE=$(head -1 $NAME) # the code of the POSCAR file
    # find the source POSCAR file in the source directory
    INFO=$(grep $CODE $DIR/3RD.POSCAR.* | head -1)
    [ x"$INFO" == x ] && echo "NOT FOUND for job-$NUM" && continue
    SOURCE_NUM="${INFO%:*}" # the number of the source POSCAR file
    SOURCE_NUM="${SOURCE_NUM##*.}"

    echo "$DIR/job-$SOURCE_NUM/ --> job-$NUM/"
    # diff 3RD.POSCAR.$NUM $DIR/3RD.POSCAR.$SOURCE_NUM
    rsync -a $DIR/job-$SOURCE_NUM/vasprun.xml ./job-$NUM/
    done
  • 针对大量POSCAR批量建立文件夹

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    #module load spglib/master
    #thirdorder_vasp.py sow 3 3 1 -10
    #生成2*2*2的超胞 0.6表示截断半径为6埃。
    #sow a b c m abc表示生成的超胞大小,m为正数的时候表示截断半径(单位:nm)(例如0.6,表示截断半径为6.0埃);为负数时表示第几近邻截断(此时应取负整数,例如-3表示第三近邻)

    for i in 3RD.POSCAR.*
    do
    s=$(echo $i| cut -d "." -f 3)
    d=job-$s
    mkdir $d
    cp $i $d/POSCAR
    cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d
    pushd $d
    popd
    done
  • 批量提交任务

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    #!/bin/sh
    root_path='pwd'

    #在每个job文件夹里进行静态计算。
    for i in {001..1050} #按顺序从job-0011050做静态vasp计算,根据生成的文件夹设定
    do
    cd job-$i
    echo "job-$i", time
    sbatch sub.sh
    wait
    cd ..
    done

    #上述静态计算完成后提取三阶力常数文件
    #cd ${root_path}

    #find job* -name vasprun.xml| sort -n| python thirdorder_vasp.py reap 2 2 1 -12 #最后四个数字跟你扩胞时一致

脚本赋予权限:chmod +x .sh;运行脚本:./sh

ShengBTE获得热输运性质

运行ShengBTE计算热输运性质

  • 获得二阶力常数计算:
    计算流程与声子谱计算流程相同,用DFPT和有限位移法都可以。将得到的FORCE_CONSTANTS重命名为FORCE_CONSTANTS_2ND。
  • 获得三阶力常数计算:见上,输出为FORCE_CONSTANTS_3RD文件。
  • 准备CONTROL文件:格式如下:
    {.line-numbers}
    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
    #CONTROL for ShengBTE
    &allocations
    nelements=2 #元素种类
    natoms=6 #原子个数
    ngrid(:)=20 20 5 #计算网格密度
    &end
    &crystal
    lfactor=0.1 #软件长度单位换算
    lattvec(:,1)=3.1613337552503711 0.0000000000000025 0.0000000000000000 #原胞晶胞参数
    lattvec(:,2)=-1.5806668776251651 2.7377953418477841 0.0000000000000000
    lattvec(:,3)=-0.0000000000000000 0.0000000000000000 12.3318234981499426
    elements= "Mo" "S" #原子顺序
    types= 1 1 2 2 2 2 #原子个数,2代表第二种元素,42代表有4个第二类原子
    positions(:,1)= 0.3333333449999998 0.6666666890000030 0.2500000000000000 #原子坐标,按顺序
    positions(:,2)=0.6666666240000012 0.3333333229999980 0.7500000189999980
    positions(:,3)= 0.3333333449999998 0.6666666890000030 0.6230719397502427
    positions(:,4)=0.6666666240000012 0.3333333229999980 0.3769280602499874
    positions(:,5)=0.6666666240000012 0.3333333229999980 0.1230720067501656
    positions(:,6)=0.3333333449999998 0.6666666890000030 0.8769280222497613
    scell(:)=4 4 1 #扩胞倍数,和前面一致
    &end
    &parameters
    T=300 #温度
    &end
    &flags
    nonanalytic=.TRUE.
    nanowires=.FALSE.
    convergence=.TRUE. #迭代法,默认最大1000
    &end

详细参数可参考手册。

  • 提交任务
  • 获得不同温度的数据,只需要依次修改CONTROL文件的T tag。或者通过脚本:
{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
for i in {350..1350..100} # 对i赋值,区间为300900,每搁10取一个数,即(300 310 320 330 ....890 900),可以根据自身情况自行设置
do #开始循环
mkdir $i #建立单一任务文件夹
cp sub.sh CONTROL FORCE_CONSTANTS_2ND FORCE_CONSTANTS_3RD $i
#将二阶三阶力学常数文件和CONTROL文件复制到单一任务文件夹
#注意,这里的CONTROL文件夹是删减过的
cd $i #进入单一文件夹内
echo \&parameters >> CONTROL #补充CONTROL文件
echo T= $i >> CONTROL #补充CONTROL文件
#温度设置 如果同时设置多个温度,同一行只会算第一个并报错。如果多次设定一个温度到同一个CONTROL文件,只会运算最后一个。
#echo scalebroad=0.1 >> CONTROL #补充CONTROL文件
echo \&end >> CONTROL #补充CONTROL文件
echo \&flags >> CONTROL #补充CONTROL文件
echo nonanalytic=.TRUE. >> CONTROL #补充CONTROL文件
echo nanowires=.FALSE. >> CONTROL #补充CONTROL文件
echo convergence=.True. >> CONTROL
echo \&end >> CONTROL #补充CONTROL文件
#这里补充CONTROL文件用的是笨方法,如果有相应改进方法可以自行更改,简洁代码
sbatch sub.sh #执行ShengBTE计算,这里的执行命令应视情况而定。
cd .. #返回初始目录
done #结束任务

数据处理

其他参考:https://www.bilibili.com/read/cv22090733/

注意事项

  • 计算得到的热导率和文献报道有较大偏差,可能偏大或者偏小,可能的原因有:

    1.采用的2阶力常数FORCE_CONSTANT中超胞内原子排序不是按照ucatomsc_zsc_y*sc_x排列;力常数文件与control文件中原胞内原子次序是否一致;

    2.shengBTE的control文件中晶格参数单位默认nm,有转换系数的参数;

    3.使用QE利用phonopy计算得到的FORCE_CONSTANT文件中,单位需要转换,不能直接利用shengBTE计算,vasp利用phonopy则单位正确;

    1. 3阶力常数的计算量较大,但精度要求很高,如0.0003Rh/bohr^2量级的力的误差,也会导致最终热导率的数值出现大的偏离;如使用cp2k计算力常数,通常使用QE等平面波软件,对于超胞体系,随着扩胞,在较大的截断半径下,所需计算量非常大,使用CP2K中的QS模块也能进行第一性原理计算,具有相对较高的计算效率,但是CP2K不同于平面波的QE等,计算中仅有gamma点计算,需要采用实空间较大超胞,如果使用和QE等相同的超胞,则会导致原子受力精度不高,一般在~0.0003Rh/bohr^2量级的误差,也会对某些体系,特别是二维体系的热导率有很大的影响。在使用cp2k进行3阶力常数大量计算之前,先进行超胞收敛测试,可以在相同的位移构型下,不断扩胞,判断力计算的变化情况和误差。用QE等计算也建议先进行k点等计算参数的收敛判断,否则会浪费大量时间,热导率计算结果却有很大偏离。

    2. 使用经验势能模型,在通过lammps进行热输运计算之前,可以先利用分子静力学软件GULP结合shengBTE模块,预判势能应用于热输运计算的准确性。但经验势能用于热输运精度和前景一般,机器势能表现更好。

    原文链接:https://blog.csdn.net/odin_linux/article/details/130084321

  • 静态计算可以各个job文件夹同时分别计算,不用按顺序来。等所有job文件夹计算完后,自行运行“find job…” 提取三阶力常数的命令。严格来说,胞的大小和三阶力常数计算的截断半径都应该进行测试,测试以所计算得到的热导率的收敛进行判断。计算资源有限无法测试的话,尽可能取大一点。Born有效电荷和介电常数的计算应该不是必要的,如果你要考虑在Gamma点附近LO-TO劈裂的话(主要针对离子晶体),则需要进行计算。可以加上,感觉要好点。
    二阶力常数和三阶力常数所用的超胞大小可以不统一,我看文献中有这种情形,不过一般作者采用的二阶力常数的超胞都是大于等于三阶力常数超胞的大小。
    原文链接:https://zhuanlan.zhihu.com/p/147810223

  • 官方网站:https://www.shengbte.org/

四声子过程

除了ShengBTE的CONTROL文件中的常规输入,FourPhonon还需要四阶力常数和一些新的标签在CONTROL文件中。请检查ShengBTE网站以获取其他标签的定义。

并行环境

1.1和1.0版本是为MPI并行性编写的。从1.2版本开始,支持四声子散射的迭代求解器,我们已经迁移到OpenMP,以处理该迭代求解器所需的大内存。未来我们可能会支持MPI+OpenMP混合并行性,以允许更多灵活性。确保在编译时添加-qopenmp,或者:

export FFLAGS=-qopenmp -traceback -debug -O2 -static_intel

四阶IFCs文件:FORCE_CONSTANTS_4TH

该文件包含四阶原子间力常数矩阵,并使用稀疏描述来节省空间。要构建四阶力常数,可以参考Fourthorder的python脚本。该力常数的格式是对三阶力常数的直接扩展,包含nb个这样的矩阵块:

  • 一行空白
  • 一个从1开始的顺序索引(从1到nb)
  • 一行表示第二个晶胞的笛卡尔坐标(单位:Å)
  • 一行表示第三个晶胞的笛卡尔坐标(单位:Å)
  • 一行表示第四个晶胞的笛卡尔坐标(单位:Å)
  • 一行表示四个相关原子的从1开始的索引,每个原子索引从1到natoms
  • 81行的力常数矩阵(单位:eV Å⁴)。开头的索引用于标记笛卡尔坐标轴。
    该文件的一个示例:
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

1
0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
1 1 1 1
1 1 1 1 -42.0002584218
1 1 1 2 0.0000000000
1 1 1 3 0.0000000000
1 1 2 1 0.0000000000
1 1 2 2 40.3875689003
1 1 2 3 0.0000000000
1 1 3 1 0.0000000000
1 1 3 2 0.0000000000
1 1 3 3 40.3875689003
1 2 1 1 0.0000000000
1 2 1 2 40.3875689003
1 2 1 3 0.0000000000
1 2 2 1 40.3875689003
1 2 2 2 0.0000000000
1 2 2 3 0.0000000000
1 2 3 1 0.0000000000
1 2 3 2 0.0000000000
1 2 3 3 0.0000000000
1 3 1 1 0.0000000000
1 3 1 2 0.0000000000
1 3 1 3 40.3875689003
1 3 2 1 0.0000000000
1 3 2 2 0.0000000000
1 3 2 3 0.0000000000
1 3 3 1 40.3875689003
1 3 3 2 0.0000000000
1 3 3 3 0.0000000000
2 1 1 1 0.0000000000
2 1 1 2 40.3875689003
2 1 1 3 0.0000000000
2 1 2 1 40.3875689003
2 1 2 2 0.0000000000
2 1 2 3 0.0000000000
2 1 3 1 0.0000000000
2 1 3 2 0.0000000000
2 1 3 3 0.0000000000
2 2 1 1 40.3875689003
2 2 1 2 0.0000000000
2 2 1 3 0.0000000000
2 2 2 1 0.0000000000
2 2 2 2 -42.0002584218
2 2 2 3 0.0000000000
2 2 3 1 0.0000000000
2 2 3 2 0.0000000000
2 2 3 3 40.3875689003
2 3 1 1 0.0000000000
2 3 1 2 0.0000000000
2 3 1 3 0.0000000000
2 3 2 1 0.0000000000
2 3 2 2 0.0000000000
2 3 2 3 40.3875689003
2 3 3 1 0.0000000000
2 3 3 2 40.3875689003
2 3 3 3 0.0000000000
3 1 1 1 0.0000000000
3 1 1 2 0.0000000000
3 1 1 3 40.3875689003
3 1 2 1 0.0000000000
3 1 2 2 0.0000000000
3 1 2 3 0.0000000000
3 1 3 1 40.3875689003
3 1 3 2 0.0000000000
3 1 3 3 0.0000000000
3 2 1 1 0.0000000000
3 2 1 2 0.0000000000
3 2 1 3 0.0000000000
3 2 2 1 0.0000000000
3 2 2 2 0.0000000000
3 2 2 3 40.3875689003
3 2 3 1 0.0000000000
3 2 3 2 40.3875689003
3 2 3 3 0.0000000000
3 3 1 1 40.3875689003
3 3 1 2 0.0000000000
3 3 1 3 0.0000000000
3 3 2 1 0.0000000000
3 3 2 2 40.3875689003
3 3 2 3 0.0000000000
3 3 3 1 0.0000000000
3 3 3 2 0.0000000000
3 3 3 3 -42.0002584218

CONTROL 文件设置

此 CONTROL 文件包含所有用户指定的设置和参数,包括晶体结构信息、展宽因子、q 网格、温度、函数等。

一般用法

要启用 FourPhonon 功能,您需要添加一个新的 &flags 名称列表:

1
2
four_phonon (逻辑,默认为 .false.):计算四声子相空间和四声子散射率。
four_phonon_iteration (逻辑,默认为 .false.):精确计算四声子散射并求解玻尔兹曼输运方程(BTE)。这仅在 four_phonon=.true. 时有效。

标志的实际用法
以下是这些标志的一些实用组合:

1
2
3
4
onlyharmonic=.true. 和 four_phonon=.true.:仅计算四声子相空间。
convergence=.false. 和 four_phonon=.true.:在 RTA 层面计算三声子和四声子散射的热导率;此设置将输出不同四声子通道的贡献。
four_phonon=.true. (收敛默认值为 .true.):使用三声子迭代方案计算热导率,但将四声子散射处理为 RTA 层面。
four_phonon=.true. 和 four_phonon_iteration=.true.:利用 BTE 的精确解计算三声子和四声子散射的热导率。

注意: CONTROL 文件中的所有其他参数,例如温度或 q 网格,适用于三声子和四声子过程。FourPhonon 包不支持纳米线功能。

加速 RTA 解的采样方法

为了从散射过程样本中估计散射率,请引用 [Z. Guo 等人,npj Comput. Mater. 10, 31 (2024).] 并使用相关标签。

参数设置

num_sample_process_3ph 和 num_sample_process_3ph_phase_space

num_sample_process_3phnum_sample_process_3ph_phase_space (整数,默认值为 -1):表示从每个模式中提取的样本数,以用于估计三声子相空间和散射率。 -1 表示不进行取样,执行严格计算。注意,当 convergence=.true.,即使用迭代方案计算三声子散射时,num_sample_process_3ph 必须为 -1,因为取样方法基于弛豫时间近似。

num_sample_process_4ph 和 num_sample_process_4ph_phase_space

num_sample_process_4phnum_sample_process_4ph_phase_space (整数,默认值为 -1):表示从每个模式中提取的样本数,以用于估计四声子相空间和散射率。 -1 表示不进行取样,执行严格计算。注意,当 four_phonon=.false. 时,num_sample_process_4phnum_sample_process_4ph_phase_space 必须均为 -1。

示例

计算三声子散射率(精确计算,迭代),并在 RTA 层面采用 100,000 个样本过程计算四声子散射:

1
2
3
4
5
6
7
8
&parameters
T=300
scalebroad=1
num_sample_process_3ph_phase_space = -1
num_sample_process_3ph = -1
num_sample_process_4ph_phase_space = 100000
num_sample_process_4ph = 100000
&end

TDEP IFCs 的输入

对于 TDEP 格式的谐波力常数,请使用以下描述的标签/脚本:

tdep=.true. 启用 TDEP 格式谐波 IFCs 的处理,该输入文件名应始终为 infile.forceconstant。如果系统是极性材料,除了常规参数(如 Born 有效电荷和介电常数)外,我们还需要用户手动输入耦合参数 (\lambda_{Ewald}),该参数在 TDEP IFCs 的末尾近似写出,如下所示:

1
0.675614929199219       # Coupling parameter in Ewald summation

注意:输入的谐波 IFCs 是未带极性修正的裸 IFCs,以便 FourPhonon 能够处理并在内部添加修正。该模块在很大程度上依赖 TDEP 源代码,我们无法保证成功。建议用户检查处理后的频率和三声子散射率。

对于非极性材料,我们还开发了一种脚本,将 TDEP 格式转换为 Phonopy 格式,名为 tdep2phonopy.f90。这样您无需激活 tdep 标志。编译 Fortran 后,运行可执行文件(以下是示例):

1
./tdep2phonopy outfile.forceconstant 10 10 1 2

这将 TDEP IFCs (outfile.forceconstant)转换为 Phonopy 格式,构建一个来自双原子单位胞的 (10 \times 10 \times 1) 超胞。请检查转换后的声子色散。我们建议使用足够大的超胞来存储所有的第二阶 IFC 中的成对相互作用。

示例

计算三声子散射率(精确计算,迭代),并在 RTA 层面计算四声子散射。使用的 TDEP 格式谐波 IFCs 并应用非解析修正,Ewald 参数为 (\lambda =0.675614929199219)。同样需要 Born 有效电荷(此处未列出),与 ShengBTE 对极性系统的处理相似。

1
2
3
4
5
6
7
8
9
10
11
&parameters
T=300
scalebroad=0.1
Ewald=0.675614929199219
&end
&flags
tdep=.TRUE.
nonanalytic=.TRUE. ! 此标签现在需要明确指定
convergence=.TRUE.
four_phonon=.TRUE.
&end

TDEP 格式的 3阶/4阶 IFCs

对于 TDEP 格式的 3阶/4阶 IFCs,我们开发了一个独立的脚本,将其转换为 ShengBTE 格式,名为 tdep2ShengBTE.f90。我们需要在同一目录下有 infile.ucposcar。编译 Fortran 后,运行可执行文件:

示例

1
./tdep2ShengBTE 3 2 outfile.forceconstant_thirdorder

这将 TDEP 的 3 阶 IFCs(outfile.forceconstant_thirdorder)转换为 ShengBTE 格式。

Output files

Besides the routine output files from previous ShengBTE program, FourPhonon generates these output files:

  • BTE.Numprocess_4ph: number of allowed four-phonon scattering processes, for each irreducible q point and phonon band
  • BTE.P4: phase space available for four-phonon processes, for each irreducible q point and phonon band
  • BTE.P4_total: total volume in phase space available for four-phonon processes
  • BTE.P4_plusplus*, BTE.P4_plusminus*, BTE.P4_minusminus*: similar to BTE.P4 but only includes contributions from ++/+-/- - processes

*Note for four-phonon scatterings, there are three different channels: recombination (++), redistribution (+-) and splitting (- -) processes.

Under temperature-dependent directories: (the unit of scattering rates is ps−1 if not specified)

  • BTE.WP4: weighted phase space available for four-phonon processes ([rad/ps]$^{-5}$, 2nd column) vs angular frequency (rad/ps, 1st column) for each irreducible q point and phonon band
  • BTE.WP4_plusplus*, BTE.WP4_plusminus*, BTE.WP4_minusminus*: similar to BTE.WP4 but only includes contributions from ++/+-/- - processes
  • BTE.w_3ph: three-phonon scattering rates for each irreducible q point and phonon band, this file replaces the original output BTE.w_anharmonic. Absorption and emission processes are written out into BTE.w_3ph_plus and BTE.w_3ph_minus
  • BTE.w_4ph: four-phonon scattering rates for each irreducible q point and phonon band. Similarly, we provide the contributions from different channels: BTE.w_4ph_plusplus, BTE.w_4ph_plusminus and BTE.w_4ph_minusminus
  • BTE.w_4ph_normal: four-phonon scattering rates of normal processes, for each irreducible q point and phonon band. This file has 5 columns and each one represents: angular frequency in rad/ps, recombination channel, redistribution channel, splitting channel, overall scattering rates from normal processes
  • BTE.w_4ph_Umklapp: four-phonon scattering rates of Umklapp processes, for each irreducible q point and phonon band. The format is the same as BTE.w_4ph_normal
  • BTE.kappa*: thermal conductivity results are written out as usual

声子谱

声子谱的计算主要有两种方法,一种是直接法,另一种是微扰密度泛函方法(DFPT)。

  • 直接法,或称frozen-phonon方法,是通过在优化后的平衡结构中引入原子位移,计算作用在原子上的Hellmann-Feynman力,进而由动力学矩阵算出声子色散曲线。直接法的缺陷在于它要求声子波矢与原胞边界supersize
    正交,或者原胞足够大使得Hellmann-Feynman力在原胞外可以忽略不计。这使得对于复杂系统,如对称性高的晶体、合金、超晶格等材料需要采用超原胞。超原胞的采用使计算量急剧增加,极大的限制了该方法的使用。
  • 1987年,Baroni、Giannozzi和Testa提出了DFPT方法。DFPT通过计算系统能量对外场微扰的响应来求出晶格动力学性质。该方法最大的优势在于它不限定微扰的波矢与原胞边界正交,不需要超原胞也可以对任意波矢求解。Castep、Vasp等采用的是一种linear response theory 的方法(或者称为 density perturbation functional theory,DFPT),直接计算出原子的移动而导致的势场变化,再进一步构造出动力学矩阵。进而计算出声子谱。

DFPT计算步骤

  • 结构优化
    声子谱的计算需要对原胞结构做高精度的充分优化,否则很容易出现虚频;
    {.line-numbers}
    1
    2
    EDIFF  = 1E-08
    EDIFFG = -1E-06
    可以采用分步优化,直到达到精度。然后优化后的结果,cp CONTCAR POSCAR:

$ phonopy -d --dim="4 4 1"

之后会生成一堆文件,其中有 SPOSCAR 和POSCAR-00* ,分别对应不同的方法密度泛函微扰法 (DFPT) 以及有限位移法。

{.line-numbers}
1
2
3
4
5
6
$ ls
phonopy_disp.yaml POSCAR-002 POSCAR-005 POSCAR-008 POSCAR-011 POSCAR-014 SPOSCAR
POSCAR POSCAR-003 POSCAR-006 POSCAR-009 POSCAR-012 POSCAR-015
POSCAR-001 POSCAR-004 POSCAR-007 POSCAR-010 POSCAR-013 POSCAR-016
$ cp POSCAR POSCAR-unitcell
$ cp SPOSCAR POSCAR

KPOINTS需适当减小,可以的话最好再进行一次收敛测试
INCAR
{.line-numbers}
1
2
3
4
5
6
7
8
9
10
IBRION =  8            (determines the Hessian matrix using DFPT)
EDIFF = 1E-08 (SCF energy convergence; in eV)
PREC = Accurate (Precision level)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
IALGO = 38 (Davidson block iteration scheme)
LREAL = .FALSE. (Projection operators: false)
LWAVE = .FLASE. (Write WAVECAR or not)
LCHARG = .FLASE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)
NSW = 1

上面参数 IBRION = 8 NSW = 1是注意修改的。
提交VASP计算

后处理
phonopy --fc vasprun.xml

读取vasprun.xml,生成FORCE_CONSTANTS文件。

编辑band.conf文件,该文件给出了高对称点路径的信息
Note:

{.line-numbers}
1
2
3
4
5
ATOM_NAME = * *  #元素名称,顺序与POSCAR保持一致
DIM = 4 4 1 #建立超胞的大小,--dim后面的
BAND = 0.5 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.5 0.0 0.0 #高对称点路径,每组高对称点之间用两个空格隔开
FORCE_CONSTANTS = READ #表示读取力常数文件FORCE_CONSTANTS
FC_SYMMETRY = .TRUE.

准备mesh.conf文件,如下所示:
{.line-numbers}
1
2
3
ATOM_NAME = Si
DIM = 2 2 2
MP = 24 24 24

{.line-numbers}
1
2
3
4
phonopy --dim="4 4 1" -c POSCAR-unitcell -p -s band.conf
(单位是THz)或者
phonopy --dim="4 4 1" --factor=524.147 -c POSCAR-unitcell -p -s band.conf
(单位是cm-1)

(读取band.conf、FORCE_CONSTANTS、POSCAR-unitcell,输出band.yaml、phonopy.yaml以及pdf图)

1
2
3
phonopy-bandplot  --gnuplot> phono.dat
如果上述命令不行,可以用旧版本phonopy
bandplot --gnuplot> phonon.dat

读取band.yaml,就可以获得二阶力常数矩阵声子谱的数据,其中第一行为高对称点的坐标。

准备mesh.conf以计算态密度和PDOS:

{.line-numbers}
1
2
3
4
5
6
ATOM_NAME = 
DIM = 4 4 1
MP=12 12 1
PDOS=1 2, 3 4, 5 6
FORCE_CONSTANTS = READ
BAND_POINTS = 71

运行命令:phonopy –dim='4 4 1' -c POSCAR-unitcell mesh.conf或者乘上因子修改单位。phonopy –dim='4 4 1' --factor=524.147 -c POSCAR-unitcell mesh.conf
得到projected_dos.dat文件,拖到Origin中画图,得到声子态密度图像(从第二列数据开始,每一列对应每个原子的态密度),对比着之前计算得到的声子谱可以分析声学支与光学支之间的相互作用关系

有限位移方法

这里需要准备和 POSCAR-00* 数目一样的文件夹,并把 POSCAR-00* 复制到每个文件夹中并命名为 POSCAR 。每个文件夹内都需要准备 POTCAR 、 POSCAR 、 INCAR 、 KPOINTS 。并且分别计算。
可用以下sh脚本快速进行:手动修改第一行文件夹数目和vasp运行文件名称。

{.line-numbers}
1
2
3
4
5
6
7
8
9
for i in {001..006}
do
mkdir $i
cd $i
cp ../POSCAR-$i POSCAR
cp ../{INCAR,KPOINTS,POTCAR,sub.sh} .
运行提交任务命令
cd ../
done

INCAR 可以设置为如下内容:
INCAR设置如下(静态计算):

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
PREC = Accurate
IBRION = -1
NSW = 1
ENCUT = 500
EDIFF = 1.0e-08
EDIFFG = -0.001
ISMEAR = 0
SIGMA = 0.05
IALGO = 38
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.

里面IBRION = -1 NSW = 1注意修改

KPOINTS需适当减小,可以的话最好再进行一次收敛测试
提交VASP计算

后处理

{.line-numbers}
1
phonopy -f ./disp-*/vasprun.xml

准备band.conf

{.line-numbers}
1
2
3
4
5
ATOM_NAME = In Se
DIM = 4 4 1
BAND = 0.5 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.5 0.0 0.0
FORCE_CONSTANTS = WRITE #注意这里是 WRITE
FC_SYMMETRY = .TRUE.
{.line-numbers}
1
2
3
4
5
6
7
8
phonopy -c POSCAR-unitcell band.conf -p -s
phonopy --factor=524.147 -c POSCAR-unitcell -p -s band.conf
#旧版本phonopy
bandplot --gnuplot> phonon.dat
就可以获得二阶力常数矩阵 FORCE_CONSTANTS 和声子谱的数据
#新版本phonopy
phonopy-bandplot --gnuplot > phonon.dat

Mind

  • 如果计算完的声子谱其他地方都好,就是在Γ点有一点点虚频,那么这个材料很可能是稳定的,只是你优化做得不够好,进一步提高优化的精度可消除这一点点虚频。
  • 对于二维材料,如果在Γ点出现很小的虚频,基本可以认为这个材料是稳定的,大部分二维材料都会有此现象;尤其是VASP结合phonopy code计算二维材料的声子谱在Γ点更是容易出现虚频;使用Quantum-ESPRESSO的PWSCF和PH模块计算声子谱对于内存的需求较小一些,且对于二维材料的声子谱计算更友好一些,后期会详细介绍Quantum-ESPRESSO计算声子谱的具体步骤。
  • 声子谱是表示组成材料原子的集体振动模式。如果材料的原胞包含n个原子,那么声子谱总共有3n支,其中有3条声学支,3n-3条光学支。声学支表示原胞的整体振动,光学支表示原胞内原子间的相对振动。

计算出的声子谱有虚频,往往表示该材料不稳定。 有些情况下,我们可以利用虚频信息使不稳定的材料变得稳定。如果声子谱的一条声学支存在虚频,主要位于Γ点和M点1/2处(对应倒格矢的1/4位置)。倒格矢的1/4,对应晶格长度的4倍。我们可能需要将原胞沿上述倒格矢方向扩大四倍,进一步优化原子位置,才可能得到比较稳定的晶胞。

消虚频

  1. DFPT方法与有限位移法互换尝试。或者更换标准胞扩胞计算。
  2. 扩胞,在虚频很小时可以尝试扩大扩胞倍数。
  3. 提高精度优化后计算声子谱。声子谱对力的精度要求特别高,所以提高精度多次优化后再计算声子谱也是可行的。
  4. 修改POTCAR,尝试不同的赝势方法。(不太建议,因为需改之后所有计算都需要统一成稳定构型的赝势。不过有效)
  5. 修改POSCAR原子坐标,再优化,再计算声子谱。
    修改方法:对于DFPT计算的声子谱,使用grep cm-1 OUTCAR查看频率,寻找最大的虚频频率。然后在OUTCAR中从最后往前找到最大的虚频频率对应的坐标。
    dx dy dz对应虚频的位移,把这些坐标与POSCAR中的坐标相加或者0.1-0.9之间的倍数,形成新的POSCAR。进行优化,再计算声子谱。
  6. 不排除真的算不出来没虚频的,据我所知,SnTe就没见到谁大大方方发文章谈声子谱的(算出来虚频消除不了不好意思讲)。
  7. INCAR参数对虚频的影响:ADDGRID影响很小,LREAL有影响。截断能和K点密度在收敛情况下影响小。

结构优化消虚频

X, Y, Z为体系中原子的坐标,也就是POSCAR中的内容。 dx, dy,dz为虚频对应的原子振动方向上的位移。我们将(dx, dy,dz)这个位移乘以一个介于0-1之间的校正因子,然后跟POSCAR中的坐标加起来即可。以z方向为例,校正因子设为0.1, 微调后的z方向坐标为: 10.709980 + (-0.676634)0.1。
脚本vib_correct.py如下:https://www.bigbrosci.com/2022/04/29/A29/

{.line-numbers}
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
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
from ase.io import read, write
import os

#获取虚频开始行
l_position = 0 #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
lines = f_in.readlines()
wave_num = 0.0
for num, line in enumerate(lines):
if 'f/i' in line:
wave_tem = float(line.rstrip().split()[6])
if wave_tem > wave_num: #获取最大的虚频
wave_num = wave_tem
l_position = num+2
# ASEPOSCAR
model = read('POSCAR')
model_positions = model.get_positions()
num_atoms = len(model)
#print(model_positions)

# 获取虚频对应的OUTCAR部分
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 72227248
#print(vib_lines)
vib_dis = []
for line in vib_lines:
#model_positions = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
vib_dis.append(vib_infor)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。
model.positions = new_positions

###保存结构
write('POSCAR_new', model, vasp5=True) # POSCAR_new是微调后的结构,用于下一步的计算(别忘了把POSCAR_new改成POSCAR)。

过渡态虚频

如果过渡态计算的时候出现好几个虚频,自己直接写个python脚本也能实现。不过使用ASE可以很方便地获取POSCAR中原子数目,保存原子固定的信息。参考如下。

{.line-numbers}
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
import os

#获取虚频开始行
l_position = 0 #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
lines = f_in.readlines()
wave_num = 0.0
for num, line in enumerate(lines):
if 'f/i' in line:
wave_tem = float(line.rstrip().split()[6])
if wave_tem > wave_num: #获取最大的虚频
wave_num = wave_tem
l_position = num+2
# 读POSCAR
with open('POSCAR') as f_pos:
lines_pos = f_pos.readlines()

# 获取虚频对应的OUTCAR部分
num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()])
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 72227248

model_positions = []
vib_dis = []
for line in vib_lines:
position = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
model_positions.append(position)
vib_dis.append(vib_infor)

# 微调结构
model_positions = np.array(model_positions)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。

###保存结构
f_out = open('POSCAR_new','w')
f_out.writelines(lines_pos[:8])
f_out.write('Cartesian\n')
for i in new_positions:
f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n')
f_out.close()

phonopy声子群速度

操作步骤

  1. 首先对结构高精度优化后扩胞,计算得到没有虚频的声子谱;
  2. 对band.conf文件进行修改,并添加群速度计算相关参数GROUP_VELOCITY = .TRUE.GV_DELTA_Q = 0.01,程序需要读取本征值和本征矢,所以必须设置 EIGENVECTORS =.TRUE.,此外为了防止vaspkit在导出数据时报错,需要删除BAND_LABELS高对称路径名称,修改后band.conf文件内容如下;

而且具体使用过程中,关于执行命令还有执行文件的tag内容还有一些可替换的关系,–gv, –group_velocity 可在执行计算命令时替代执行文件中的GROUP_VELOCITY = .TRUE.内容, –gv_delta_q 可在执行计算命令时替代执行文件中的(GV_DELTA_Q)内容

  1. 对声子谱重新进行计算,打开phonopy.yaml文件可以看到,phonopy在计算声子谱的过程中考虑到了群速度;打开band.yaml文件中也保存了有关群速度的相关信息;同时在计算声子谱中phonopy所储存的数据文件mesh.yaml和band.yaml文件中也保存了有关群速度的相关信息。
  2. 随后进行数据提取,以vaspkit1.12版本为例,打开vaspkit,输入73,打开VASP2other Interface 子功能,选择739 Sort Phonon Band Structure for Phononpy,运行完成即可获得group velocity相关数据以及各原子的声子谱数据;在运行过程中我们可以看到软件说明find group velocity并写入输出文件中。
    1.5版本为789功能。

格林艾森常数(Grüneisen parameters)

为表示声子非简谐效应的参数。phonopy计算模式格林艾森参数方式见官网。
以NaCl为例,首先需要进行三个声子谱计算,分别是平衡体积(orig),稍小体积(minus)和稍大体积(plus)的结构,计算时改变晶格参数使得体积获得改变。

计算完成后获得不同体积下的原子间力常数文件FORCE_SETS 或 FORCE_CONSTANTS (后续命令中需要添加–readfc参数)。

随后通过phonopy-gruneisen(旧版gruneisen)命令获得对应于声子带结构中倒易空间中的路径的模式格林艾森常数(具体设定可根据自身计算情况更改)
phonopy-gruneisen orig plus minus --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --band="1/2 1/4 3/4 0 0 0 1/2 1/2 1/2 1/2 0.0 1/2" -p -c POSCAR-unitcell
在倒易网格中的模式格林艾森常数(具体设定可根据自身计算情况更改):
phonopy-gruneisen orig plus minus --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --mesh="20 20 20" -p -c POSCAR-unitcell --color="RB"
phonopy-gruneisen 命令选项包括–dim,–mp, –mesh,–band,–pa, –primitive_axis,–readfc,–band_points,–nac,–factor,–nomeshsym,-p,-c,-s, –save,-o等,功能上均与phonopy命令相同,但除上述列出的选项外,可以在phonopy命令上使用的其它选项是不能在phonopy-gruneisen 上使用的。如有必要,可以自行修改phonopy-gruneisen代码。

如查看帮助,可输入
phonopy-gruneisen -h
通过 phonopy-gruneisenplot 命令导出计算得到的数据,使用方式类似于phonopy-bandplot

原文链接:https://blog.csdn.net/icehoqion/article/details/130550807

热学性质

参考:
http://supersunsir.com/2020/08/01/vasp/vasp_phonopy/

热膨胀系数

我们通过准简谐近似(QHA)的方法计算热膨胀系数和等压热容等性质,

  • 首先,我们要把优化好的POSCAR加应变(比如0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04的缩放系数,这样就有9个结构文件),
  • 然后 把每一个POSCAR都去重复步骤(2)也就是去计算每一个缩放系数声子谱。大概就是这个效果
  • 最后我们通过phonopy里面的QHA模块 去得到此材料的非偕性质phonopy-qha -p -s e-v.dat *00/thermal_properties.yaml就可以得到很多热力学性质比如热膨胀系数,格林爱森参数,自由能,等压热容等等。
  1. 先进行高精度的晶格优化。
  2. 修改缩放系数为0.95-1.05新建10个POSCAR.
  3. 对11个POSCAR晶体进行声子谱计算
  4. thermo.conf文件后处理生成11个对应的thermal_properties.yaml
  5. QHA计算热膨胀系数

另一个写的步骤:

  1. 高精度的原胞优化
  2. 在此基础上施加形变(梯度可以设置为千分之一)
  3. 优化多个结构
  4. 使用DFPT方法计算每个结构的声子谱(如果结构有问题,使用phonopy –symmetry)
  5. 在每个文件夹中得到thermal_properties.yaml
  6. 运行脚本得到e-v.dat
  7. 运行phonopy-qha的相关命令得到各种dat文件
  8. Over

可以参考作者的thermo.conf文件:

{.line-numbers}
1
2
3
4
5
6
7
DIM = 2 2 2
ATOM_NAME = sio
MP = 31 31 31
TPROP = .TRUE.
TMAX = 2000
TSTEP = 2
FORCE_CONSTANTS = READ

然后是如何写e-v.dat文件

  • 大概意思就是找到之前计算的那11个POSCAR的vasprun.xml结果文件,然后记事本打开搜索最后一个e_wo_entrp。
  • 下载phonopy官方文件中Si-QHA的example,
  • 解压vasprun_xmls.tar.lzma我们可以得到算例中用到的vasprun.xml文件。

不同温度下热力学性质计算(自由能、热容、振动熵)

  • 在DFPT或者Finite displacement方法计算声子谱的基础上,得到mesh.yaml文件

  • 新建一个文件夹,将FORCE_CONSTANTS、POSCAR、mesh.conf放入该文件夹中(POSCAR为原胞)

  • 得到mesh.yaml文件

    {.line-numbers}
    1
    phonopy mesh.conf -c POSCAR

    mesh.conf的内容

    {.line-numbers}
    1
    2
    3
    4
    ATOM_NAME = Hf Si O
    DIM = 2 2 2
    MP = 13 13 13
    FORCE_CONSTANTS = READ

    得到thermal_properties.yaml文件

phonopy -c POSCAR mesh.conf -t
得到thermal_properties.dat文件
phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat
将dat文件导入到Origin中画图

dat文件的5列数据分别为:

{.line-numbers}
1
2
3
4
5
temperature[K]
Free energy [kJ/mol]
Entropy [J/K/mol]
Heat capacity [J/K/mol]
Energy [kJ/mol]

下面是能用到的各种命令

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
扩胞用
phonopy -d --dim="2 2 2" -c POSCAR_UNIT
查找晶格对称性
phonopy --symmetry
收集力常数
phonopy --fc vasprun.xml
得到band.yalm文件(好像不需要这一步)
phonopy --dim="2 2 2" -c POSCAR-UNIT -p band.conf
得到thermal_properties.yaml文件
phonopy -c POSCAR mesh.conf -t

下面是产生e-v.dat文件的脚本,需要根据自己的目录名进行修改:

{.line-numbers}
1
2
3
4
5
6
7
8
9
for i in -0.005 -0.004 -0.003 -0.002 -0.001 0.001 0.002 0.003 0.004 0.005
do
cd AC_Uni$i;
pwd;
volume=`grep volume vasprun.xml | tail -1| awk '{printf("%20.13f", $3)}'`;
energy=`grep 'E0' OSZICAR | tail -1 | awk '{printf "%12.4f \t", $5}'`;
echo ${volume} ${energy} >> ../e-v.dat;
cd ../; #Exit the Strain folder
done

声子谱振动可视化

  1. 上面生成band.yaml文件,利用JMOL打开。https://jmol.sourceforge.net/
    或者网站:https://henriquemiranda.github.io/phononwebsite/phonon.html

  2. band.conf添加标签ANIME=0 0 0将000这个Q点的振动模式输出到anime.ascii,用v_sim软件打开即可。https://l_sim.gitlab.io/v_sim/index.en.html

  3. 运行命令 phonopy -f vasprun.xml-{001,002}得到FORSETS
    运行phonopy –dim=”3 3 3” –ct=”0 0 0”得到gamma点的振动模式密度

  4. DFPT方法得到的OUTCAR可以利用PyVibMS、JMOL

  5. 从band.yaml生成vesta读取:https://github.com/AdityaRoy-1996/Phonopy_VESTA

注意不要输出声子群速度,写出声子本征矢量。需要POSCAR.vesta格式。

从OUTCAR生成vesta读取:https://github.com/Stanford-MCTG/VASP-plot-modes

  1. 其他参考:https://www.bilibili.com/read/cv11848434/

声子谱反折叠

https://github.com/yuzie007/upho

热位移

热学性质是根据倒易空间采样网格上的声子频率计算得出的,必须设置Mesh tag进行使用,并且必须检查它们相对于取样网格的收敛性。 通常这种计算对计算量要求不高,因此随着采样网格密度的增加,收敛性很容易达到。此外,计算时Phonopy对频率的默认单位为Thz。
在DFPT和冷冻声子的计算基础上,在thermal_properties.conf 中设置:

1
2
3
4
5
6
TPROP = .TRUE. #进行热学性质计算
DIM = 1 1 1 #扩胞情况
MP = 16 16 8 #Mesh网格密度
TMAX = 2000
TMIN = 0
TSTEP = 10 #温度上下限默认为100,0,step默认值为10

计算得到thermal_properties.yaml文件:

phonopy -c POSCAR thermal_properties.conf -t

得到thermal_properties.dat文件:

phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat

五列数据分别为:temperature[K]、Free energy [kJ/mol]、Entropy [J/K/mol]、Heat capacity [J/K/mol]、Energy [kJ/mol]。

Phonopy官网还给出了两个参数:
PRETEND_REAL = .TRUE.
这里将虚频声子当成实声子进入计算,仅在有小虚频实作测试用;
CUTOFF_FREQUENCY = 0.1
默认值为0,将虚频声子频率视为负值,计算高于截止频率的声子。

热位移计算

计算投影到笛卡尔轴上的均方位移与温度的函数关系。声子频率以 THz 为单位(phononopy 的默认设置),用于获得均方位移。同样需要与Mesh tag共同使用,温度设置与前述相同。新建thermal_displacements.conf用于输入tag。

1
2
3
4
5
6
TDISP = .TRUE. #进行均方位移计算
DIM = 1 1 1 #扩胞情况
MP = 16 16 8 #Mesh网格密度
TMAX = 2000
TMIN = 0
TSTEP = 10 #温度上下限默认为100,0,step默认值为10

XYZ_PROJECTION = .TRUE.
打开该tag,将对特征向量按实空间进行投影,得到不同原子在a,b,c三个轴向上的位移投影,数据存储为frequency atom1_x atom1_y atom1_z atom2_x atom2_y atom2_z…形式。

PROJECTION_DIRECTION = 1 1 0
打开该tag,将原子位移在指定实空间轴向上进行投影。
phonopy -c POSCAR thermal_displacements.conf -p -s
执行后,输出结果将保存在thermal_displacements.yaml中。注意,计算的原子位移为均方位移。

此外的两个tag:

TDISPMAT = .TRUE.
打开该tag时,将以3*3矩阵形式储存原子在六个方向上的均方位移xx, yy, zz, yz, xz, xy,结果写入thermal_displacement_matrices.yaml中。

TDISPMAT_CIF = 1273.0
该tag用于指定计算在某温度(1273K)下,原子的均方位移矩阵,结果输出在tdispmat.cif中。

Phonopy其他功能

  • phonopy --symmetry -c POSCAR-unit

BPOSCAR和PPOSCAR生成,分别是标准胞和原胞,vasp格式。

  • phonopy --rd 10 --dim 2 2 2 --pa auto --amplitude 0.03 -c POSCAR-unitcell

–rd 10 代表生成10个随机位移结构,–dim 2 2 2扩胞倍数,–pa auto沿三个坐标轴投影方向,–amplitude 0.03位移幅度,-c指定对象。

输出”phonopy_disp.yaml” and supercells

  • phonopy-load phonopy_disp_orig.yaml --rd 1000 --temperature 300

generate 1000 supercells with displacements at 300K.

  • New in v2.25 Symfc (symfc/symfc) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use symfc, its python module has to be installed via conda-forge or pip.

New in v2.3 ALM (ttadano/ALM) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use ALM, its python module has to be installed via conda-forge or building it. The installation instruction is found here.

When ALM is used, please cite the paper: T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018).

FC_CALCULATOR = ALM

参考:

https://zhuanlan.zhihu.com/p/355317202

https://mp.weixin.qq.com/s?__biz=MzI3MzIyMjI2Nw==&mid=2247488523&idx=1&sn=90a544712f166e248993e1198eb40981&chksm=eb27cb1fdc50420921bb11d3bc8f4833c14418b0d0c0d755af6687fca024fc9d15d999ae7813&scene=21#wechat_redirect

https://blog.sciencenet.cn/blog-2909108-1154273.html

https://zhuanlan.zhihu.com/p/642116307

https://zhuanlan.zhihu.com/p/25432728

一、安装编译

wannier90安装

下载: https://wannier.org/download/

2.1版本是可以与VASP结合的最新版本。但是2.1版本默认安装与VASP接口并不好,主要是借助肖承诚博士写了一个Fortran的接口(https://github.com/Chengcheng-Xiao/VASP2WAN90_v2_fix), 需要注意的是这个接口是针对VASP 5.4.4版本的。1.2版本与VASP的接口是好的,所以1.2版本默认安装就好。对比两个版本的不同,1.2版本的主要缺点就在于不能构建向上自旋和向下自旋的能带,也就是没有自旋轨道耦合作用的铁磁和反铁磁体系。 0201

1.1 wannier90 2.1安装

首先,解压安装包

tar -zxvf wannier90-2.1.0.tar.gz

其次,进入文件夹

cd wannier90-2.1.0/

然后,准备编译文件(这里老王用的是ifort,注意要检查ifort和mpiifort执行命令)

cp config/make.inc.ifort make.inc

接着,编译

make

完成后,编译库,得到libwannier.a文件

make lib

1.2 VASP编译(注意以上接口VASP2WAN90_v2_fix只针对VASP 5.4.4)

首先拷贝VASP2WAN90_v2_fix接口文件中的mlwf.patch 到VASP代码编译目录下。然后执行如下命令。

patch -p0 < mlwf.patch

接着在VASP makefile.include 文件中加入下面两行。注意路径。

1
2
3
CPP_OPTIONS+=-DVASP2WANNIER90v2

LLIBS+=/path/to/your/wannier90_distro/libwannier.a

最后编译即可

make all

wannier-tools安装

开源软件包WannierTools,是一个用于研究新型拓扑材料的软件。此代码在紧束缚模型中工作,可以由另一个软件包Wannier90生成。

它可以通过计算Wilson loop来帮助对给定材料的拓扑相进行分类,通过角分辨光电发射(ARPES)和红外光谱技术,可以得到表面态光谱扫描隧道显微镜(STM)实验。它还可以识别Weyl/Dirac点和节点线结构的位置,计算闭合动量环周围的贝里相位和部分布里渊区(BZ)的贝里曲率。
此外,WannierTools还可以计算非磁性金属和半金属的普通磁电阻利用玻尔兹曼输运理论,计算给定磁场方向和强度下的朗道能级谱;并从超胞计算中得到展开的能谱。

wanniertools官方安装教程

http://www.wanniertools.com/installation.html

安装依赖环境

1、Intel编译器即oneapi(Fortran compiler (Gfortran or ifort),MPICH version higher than 2.1.5,Lapack and Blas library, (Intel MKL recommended))

2、Arpack-ng

安装步骤

下载wanniertools和arpack安装包并分别解压

1
2
unzip wannier_tools-master.zip
tar vzxf arpack96.tar.gz

先进入arpack,修改编译文件并编译
在ARmake.inc中修改编译路径,这里以实际安装路径为准

1
home = $(HOME)/ARPACK

修改编译器,参考VASP安装时的makefile.include

1
2
FC=f77
FFLAGS= -O -cg89

改为

1
2
3
FC         = mpiifort
FCL = mpiifort -mkl=sequential
FFLAGS = -O -cg89

然后执行make编译

make lib

直至产生libarpack_SUN4.a文件,该文件会因为不同机器不同环境而不同名称,只要确定好,然后在编译wanniertools的时候指定好就行了

在wanniertools的src文件夹的makefile文件中该默认值为

1
2
3
#ARPACK LIBRARY

ARPACK=/Users/quan/work/workplace/ARPACK/libarpack MAC.a

修改为刚刚已经生成的文件路径,这里使用的是Makefile.intel-mpi,重命名为makefile

1
cp Makefile.intel-mpi Makefile

然后编译即可

1
make

最后生成wt.x文件即为安装成功

1
cp -f wt.x../bin

二、计算过程(以石墨烯为例)

1、结构优化加自洽计算

POSCAR为

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
Surf-mp-48                             
1.00000000000000
1.2344083195740001 -2.1380573715510001 0.0000000000000000
1.2344083195740001 2.1380573715510001 0.0000000000000000
0.0000000000000000 0.0000000000000000 12.0000000000000000
C
2
Direct
0.0000000000000000 0.0000000000000000 0.5000000000000000
0.3333329999917112 0.6666670000099160 0.5000000000000000

INCAR为

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
SYSTEM=Graphene
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ISPIN = 1 (Non-Spin polarised DFT)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = .FALSE. (Projection operators: automatic)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = A (Precision level)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)

Electronic Relaxation
ALGO = N
ISMEAR = 0 (Gaussian smearing; metals:1)
SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 90 (Max electronic SCF steps)
NELMIN = 6 (Min electronic SCF steps)
EDIFF = 1E-08 (SCF energy convergence; in eV)

Ionic Relaxation
NSW = 0 (Max ionic steps)
IBRION = -1 (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
ISIF = 3 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -2E-02 (Ionic convergence; eV/AA)

加SOC在自洽计算加上
LSORBIT=T

2、非自洽计算投影能带

在自洽计算基础上,cp -r scf band
修改INCAR

1
2
3
ISTART =  1         
ICHARG = 11
LORBIT = 10/11

KPOINTS 改为高对称点路径。

画一下能带结构,然后再做一下Fatband分析,了解费米面附近能带是有哪些成分构成的。以C的pz轨道。

3、非自洽vasp产生.win文件

在自洽计算基础上 cp -r scf wan1
修改INCAR

1
2
3
4
5
ISTART =  1         
ICHARG = 11
NBANDS = 40
LWANNIER90 = .TRUE.
NUM_WANN = 2

注意,INCAR中使用NPAR参数会出错

KPOINTS、POSCAR、POTCAR不改。

准备好wannier90.win文件,同时在INCAR里头添加一行LWANNIER90=.TRUE.和一行ICHARG=11,然后再运行一遍VASP。 顺利的话,会产生 wannier90.amn, wannier90.mmn, wannier90.eig, wannier90.chk 等4个文件。wannier90.win里头的内容也有可能会被改变。包括结构信息,k点信息。有了这三个文件就可以进一步调节参数,构造wannier函数了。请注意wannier90.win里头的num_bands要和OUTCAR里头的NBANDS一致。.mmn积分矩阵和.amn布洛赫态在轨道上的尝试投影和UNK0000x.1存储实空间单胞的布洛赫态相关信息,这些文件是运行链接wannier90的DFT计算输出的结果,而gaas.win是计算参数设置文件

4、wannier90.x计算

运行wannier90主程序。mpirun -np 40 wannier90.x wannier90 &

注意在第3步打开hr_plot = T,关闭write_hr = true;在第4步打开write_hr = true,关闭hr_plot = T

这里后面一个winner90代表wannier90.win的文件.win前面的部分。如果你的win文件是name.win运行的地方就可以是

wannier90.x name &

运行过程中当num_bands > num_wann 有个迭代过程,控制迭代步数可以通过

1
2
dis_num_iter  =  200   # 迭代1,默认200步
num_iter = 200 # 迭代2,默认200步

最后会生成wannier90.chk文件,restart计算会读取这个文件。

加入能带部分的话第一次运行就会出来,wannier90_band.dat文件,这个可以用来和DFT结果对比。
构造好Wannier函数以后,仔细检查一下wannier90.wout,
搜索”Final State”,在这里可以看到每条Wannier轨道的展宽,由此可以检查你构造的Wannier函数是否足够局域。如果有某个Wannier函数的展宽
比晶格常数大许多,就表明你的Wannier函数构造不是很理想,需要不断调节投影轨道,解纠缠窗口以及Frozen窗口等。从这些轨道的中心和展宽也能
看出轨道之间的简并性如何,由此也可以判断所构造的Wannier函数对称性如何,这个也会影响后续的分析,特别是表面态的计算。查看wannier90.wout末尾的Final State的最后一列小于3,其结果是可以的。

5、

紧接着,你就可以开始书写WannierTools的主输入文件了,
目前叫input.dat,以后会改成http://wt.in (WannierTools v2.2之后)。 这里不具体介绍每一个参数,在说明文档中有详细的介绍,中文版的说明文档正在翻译中。
input.dat文件有模板,使用的时候直接拷贝然后加以修改。很多内容都可以直接从wannier90.win和wannier90.wout中拷
贝。 Bi2Se3的主输入文件input.dat可以在附件中下载得到。我们可以通过控制输入文件中的一些参数来实现我们所需要的功能。

常见报错

  1. 首先跑接口会生成3个文件.amm .mmm .eig通过这三个文件构造wannier函数。也可以提前把.win文件准备好放进去。以上三个文件生成后运行wannier90.x会成.chk文件

  2. Restar=T是需要读取.chk文件的

Error: restart requested but wannier90.chkfile not found

原因:加了restart = plot 这个命令,但却没有wannier90.chk这个文件。Restart是需要读取wannier90.chk的。

解决办法:注释掉restart = plot,然后重新运行wannier90.x wannier90 &得到wannier90.chk.

  1. 问题:Wannier90 is running in LIBRARY MODE

Setting up k-point neighbours…

Ignoring in input file

解决办法:尽量串行,不要并行。请再次运行 $postw90 wannier90的为输出的数据

  1. Spinors与文件INCAR中的ISPIN的设置有关,同时存在或者同时去掉,否则影响.amn文件的形成。(自旋轨道耦合)

  2. 谨慎用use_bloch_phases=T一般是F

  3. 问题:param_read: mismatch in wannier90.eig

解决办法:每一步计算的核数要相等

  1. 问题:dis_windows: More states in the forze window than target WFs

解决办法:把dis_windows改大点

  1. Fewer projections defined than the number of wannier functions requested.

投影的轨道太少

解决办法:增加投影的轨道

  1. Kmesh_get_bvector: Not enongh bevectors found

解决办法:添加kemsh_tol

  1. Wanted band:1 found band:17

Wanted kpoint:2 found kpoint:1

解决办法:常见原因是num_bands设置错误,返回band计算文件夹查看NBANDS

  1. Dis_windows: Energy window contains fewer states than number sf target WFs

解决办法:把dis_windows改大点

  1. Dis_windows: more states in the frozen window than WFs

解决办法:调节frozen window, 投影太多或者这个窗口太大涉及到别的轨道

  1. Wannier 运行有两种模式一是library mode,一种是单机模式stand-alone mode

解决办法:library mode是需要进行串行。

  1. param_get_projection: too many projections defined

解决办法:begin projections random

原因:num_wann 与设置投影数不匹配。

解决办法:不加soc:num_wann 等于每种原子投影轨道数之和

soc:num_wann 等于两倍的每种原子投影轨道数之和

  1. Gnuplot的报错信息”arc_bulk.gnu”, line 26: warning: No unsable data in this plot to auto-scale axis range

  2. param_get_projection: Problem reading m state into string

解决办法:已解决Fe: s;px;py;pz;dz2;dxz;dyz;dx2-y2;dxy

不是Fe:s;p;d

  1. 能带拟合的不好可能是因为,投影的不对,还有核数和能带数的关系

解决办法:太马虎了,这样的问题肯定是K点写错了

能带布里渊区的宽度不一样有可能是wannier里规定的K点数并没有完全计算。

  1. num_iter = 2000 #wannier center 迭代次数

num_print_cycles = 40 #wannier center 迭代每40步输出一次

  1. dis_num_iter = 1200 #解纠缠的迭代步数#

  2. num_iter最小化过程中的总迭代次数。如果您希望生成投影wf而不是最大本地化的wf,则设置num_iter=0(请参阅教程中的示例8)?默认值为100

  3. dis_mix_ratio = 1.d0 # 0.5<dis_mix_ratio<1, 默认值0.5

#在解纠缠过程中,要利用混合参数进行收敛,值为0.5是“安全”的选择。使用1.0(即不混合)通常会得到更快的收敛,但可能会导致 最小化,在某些情况下是不稳定的。

dis_conv_tol = 1.0E-10 # Delta < 1.000E-10 over 3 iterations

  1. dis_win_max = 15

dis_froz_max = 14

num_wan = 12

exclude_bands : 1-11,26 #在使用这个参数时,num_bands = 总能带数减去排除在外的能带数

num_bands = 14

#如果使用能量窗口规定想要拟合的能带区间,需要同时使用这个几个参数

  1. wannier_plot_supercell =2 #在一个对应于“超胞”的网格上生成WFs

  2. wannier_plot_list =2 #列举要画出WFs的序号,与wannier center 的序号是一致的

  3. 运行接口的情况有三种,一是选择投影能带数发生了改变,二是,需要写入wannier90.win的输入文件;三是需要在INCAR添加关键参数,输出wannier.x需要的文件。(仅仅是调节能量窗口是不需要重新运行接口的)

  4. 无没有projection是肯定没有amn的

简介

https://github.com/hackingmaterials/Amset
https://hackingmaterials.lbl.gov/amset/
https://www.nature.com/articles/s41467-021-22440-5.pdf

劳伦斯·伯克利国家实验室Alex M. Ganose和Anubhav Jain开发了一种根据第一性原理输入来计算固态半导体和绝缘体载流子散射率的有效计算方法(AMSET)。本方法扩展了现有的极性和非极性电子-声子耦合、离子化杂质和基于各向同性能带结构的压电散射机制的公式以支持高度各向异性的材料。他们通过计算23种半导体的电子传输特性来测试这个公式,包括48个大原子$CH_3NH_3PbI_3$混合钙钛矿,并将结果与实验测量结果和更详细的散射模拟进行比较。相对于实验的斯皮尔曼迁移率秩系数(rs = 0.93)与使用恒定弛豫时间近似值(rs = 0.52)获得的结果相比有明显提高。他们发现,他们的方法可提供与最新技术类似的精度,而计算成本约为1/500,因此可在高通量计算工作流程中使用,以精确筛选载流子迁移率、寿命和热电功率。

  • Inputs obtainable from first-principles calculations. The primary input for AMSET is an uniform band structure calculation.
  • Scattering rates calculated in the Born approximation using common materials properties such as phonon frequencies and dielectric constants.
  • Transport properties calculated through the Boltzmann transport equation.
  • Efficient implementation that can run on a personal laptop.

安装

依赖

pymatgen==2023.3.23
scipy==1.10.1
monty==2022.9.9
matplotlib==3.7.1
BoltzTraP2==22.12.1
tqdm==4.65.0
tabulate==0.9.0
memory_profiler==0.61.0
spglib==2.0.2
click==8.1.3
sumo==2.3.6
h5py==3.8.0
pyFFTW==0.13.0
interpolation==2.2.4
numba==0.56.4

conda environment

pip install amset

Developer Installation

{.line-numbers}
1
2
3
git clone https://github.com/hackingmaterials/amset.git
cd amset
pip install -e .

准备工作

Structural relaxation:

“tight” calculation settings

{.line-numbers}
1
2
3
4
ADDGRID = True  
EDIFF = 1E-8
EDIFFG = -5E-4
PREC = Accurate

Dense uniform band structure and wave function coefficients:

2倍的K点密度,HSE06 is required,输出波函数,

wave ``` 得到
1
2
3
4
5
6
```wavefunction.h5 file ```
### Elastic constants:
"tight"relax
```javascript {.line-numbers}
NSW = 1
IBRION = 6

Deformation potentials:

amset deform create
single point calculation,

{.line-numbers}
1
2
NSW = 1   
ICORELEVEL = 1

amset deform read undeformed def-1 def-2 def-3
产生deformation.h5 file。

Dielectric constants, piezoelectric constants and polar-phonon frequency:

DFPT方法,Note, DFPT cannot be used with hybrid exchange-correlation functionals. In these cases the LCALCEPS flag should be used in combination with IBRION = 6

{.line-numbers}
1
2
3
NSW = 1
IBRION = 8
LEPSILON = True

amset phonon-frequency
读取DFPT的 vasprun.xml 获得 dielectric constants 和 polar phonon frequency。

运行

有两种方式,Python脚本
或者通过读取settings.yaml和输入文件加命令行运行

脚本

cp path/example/Si/Si.py .
脚本如下:

{.line-numbers}
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
import warnings
from amset.core.run import Runner
from amset.plot.rates import RatesPlotter
warnings.simplefilter("ignore")
settings = {
#general settings
"scattering_type": ["IMP", "ADP"],
"doping": [-1.99e14, -2.20e15, -1.72e16, -1.86e17, -1.46e18, -4.39e18],
"temperatures": [300],
"bandgap": 1.14,
#electronic_structure settings
"interpolation_factor": 50,
#materials properties
"deformation_potential": "deformation.h5",
"elastic_constant": [
[144, 53, 53, 0, 0, 0],
[53, 144, 53, 0, 0, 0],
[53, 53, 144, 0, 0, 0],
[0, 0, 0, 75, 0, 0],
[0, 0, 0, 0, 75, 0],
[0, 0, 0, 0, 0, 75],
[0, 0, 0, 0, 0, 75],
],
"static_dielectric": [[11.7, 0, 0], [0, 11.7, 0], [0, 0, 11.7]],
"high_frequency_dielectric": [[11.7, 0, 0], [0, 11.7, 0], [0, 0, 11.7]],
#performance settings
"write_mesh": True,
}

if __name__ == "__main__":
runner = Runner.from_vasprun("vasprun.xml", settings)
amset_data = runner.run()

plotter = RatesPlotter(amset_data)
plt = plotter.get_plot()
plt.savefig("Si_rates.png", bbox_inches="tight", dpi=400)

python Si.py 即可

命令行设置

需要准备settings.yaml文件:以Si为例:

{.line-numbers}
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
 # general settings
scattering_type: [IMP, ADP]
doping: [-1.99e14, -2.20e15, -1.72e16, -1.86e17, -1.46e18, -4.39e18]
temperatures: [300]
bandgap: 1.14

# electronic_structure settings
interpolation_factor: 50

# materials properties
deformation_potential: deformation.h5
elastic_constant:
- [144, 53, 53, 0, 0, 0]
- [ 53, 144, 53, 0, 0, 0]
- [ 53, 53, 144, 0, 0, 0]
- [ 0, 0, 0, 75, 0, 0]
- [ 0, 0, 0, 0, 75, 0]
- [ 0, 0, 0, 0, 0, 75]
static_dielectric:
- [11.7, 0, 0]
- [0, 11.7, 0]
- [0, 0, 11.7]
high_frequency_dielectric:
- [11.7, 0, 0]
- [0, 11.7, 0]
- [0, 0, 11.7]
#performance settings
write_mesh: true

amset run –directory path/to/files 运行即可
其他设置可以通过命令行加入:
amset run –interpolation-factor 20
画图:
amset plot rates mesh_99x99x99.h5

Output filesThe transport results will be written to transport_99x99x99.json. The 99x99x99 is thefinal interpolated mesh upon which scattering rates and the transport properties arecomputed.The write_mesh option is set to True, meaning scattering rates for all doping levelsand temperatures will be written to the mesh_99x99x99.h5 file.

设置相关参数:

General settings

  1. scattering_type:Default: auto ;单独设置: ADP,IMP,POP
    ADP (acoustic deformation potential scattering)
    IMP (ionized impurity scattering)
    PIE (piezoelectric scattering)
    POP (polar optical phonon scattering)
    CRT (constant relaxation time)
    MFP (mean free path scattering)
  2. doping:负值代表电子,正值代表空穴
    Default: ['1.e15', '1.e16', '1.e17', '1.e18', '1.e19', '1.e20', '-1.e15', '-1.e16', '-1.e17', '-1.e18', '-1.e19', '-1.e20']
    单独设置有两种: 1E13,1E14,1E15,1E16,1E17,1E18,1E19,1E20 或者 1E13:1E20:8
  3. temperatures:Default: [300]
    单独设置:300,400,500,600,700,800,900,1000 或者300:1000:8
  4. interpolation_factor:Default: 10
    控制插值密度,最终的k点等于能带计算的k点乘上插值因子。
  5. wavefunction_coefficients:Default: wavefunction.h5
    波函数系数读取的方式,产生该文件可以通过amset wave获得。
  6. bandgap:
    设置该参数会自动采用scissor修正,因此不能与scissor选项共用,且对金属不起作用。
  7. scissor:
    正值代表带隙打开,负值代表缩小,对金属无效。
  8. zero_weighted_kpoints:
    设置方式:
    keep: Keep zero-weighted k-points in the band structure.
    drop: Drop zero-weighted k-points, keeping only the weighted k-points.
    prefer: Drop weighted-kpoints if zero-weighted k-points are present in the calculation (useful for cheap hybrid calculations).
  9. use_projections:Default: False
    不建议设置。该参数使用投影计算波函数交叠,往往导致差的性能,使用时设置LORBIT = 11
  10. free_carrier_screening:Default: False
    控制是否屏蔽极化光学声子和压电散射率。这可以导致在高载流子浓度下散射率的大幅度降低。

Material settings

  1. high_frequency_dielectric:Required for: POP, PIE
  2. static_dielectric:Required for: IMP, POP
  3. elastic_constant:Required for: ADP, PIE
  4. deformation_potential:Required for: ADP
  5. piezoelectric_constant:Required for: PIE
  6. defect_charge:Required for: IMP。Default: 1
  7. pop_frequency:Required for: POP amset phonon-frequency获得。
  8. compensation_factor:Required for: IMP Default: 2
  9. mean_free_path:Required for: MFP
  10. constant_relaxation_time:Required for: CRT 不建议和其他散射机制设置。

Performance settings:控制速度和精度,一般默认。

  1. energy_cutoff:Default: 1.5
  2. fd_tol:Default: 0.05. 0-1之间,越小计算K点越多。
  3. dos_estep:Default: 0.01.越小越好越贵。
  4. symprec:Default: 0.01。
  5. nworkers:Default: -1。代表用所有的机器。其他设置代表几个。
  6. cache_wavefunction:Default: True。加速计算,增加内存。如果内存不够关闭。

Output settings:输出文件设置

  1. calculate_mobility:Default: True。是否计算n/p载流子迁移率。
  2. separate_scattering_mobilities:Default: True。是否计算单独的散射机制的迁移率。
  3. write_mesh:Default: False。是否将全K点结果写入磁盘。对于大的interpolation_factor可以打开。
  4. file_format:Default: json。输出文件的格式:json, yaml, and txt.对于write_mesh=True不能用TXT。
  5. write_input:Default: False。是否将输入设置写入amset_settings.yaml.
  6. print_log:Default: True。是否输出日志。

Summary of scattering rates

Mechanism Code Requires Type
Acoustic deformation potential scattering ADP n- and p-type deformation potential, elastic constant Elastic
Piezoelectric scattering PIE high-frequency dielectric constant, elastic constant, piezoelectric coefficient (e) Elastic
Polar optical phonon scattering POP polar optical phonon frequency, static and high-frequency dielectric constants Inelastic
Ionized impurity scattering IMP static dielectric constant Elastic

其他功能

you can use amset which relies on Fourier interpolation of the band
structure using BoltzTraP2. There are two options.

  • Use amset eff-mass: This calculates the conductivity effective mass as defined here: https://www.nature.com/articles/s41524-017-0013-3
  • Use amset plot band --stats vasprun.xml: This calculates the effective mass in the same way as sumo but using a band structure interpolated from a uniform calculation.

自洽计算

VASP做静态自洽的目的主要是为了得到好的电子密度和波函数文件,为了使后续的性质计算可以读取电子密度和波函数,进而增加收敛速度。
首先进行 高精度 的结构优化。然后复制文件, cp CONTCAR POSCAR 做一次自洽计算。
自洽计算要求有密的k网格点,在计算条件允许的情况下要求自洽的k网格点大致为优化时的2倍左右。
修改标签

{.line-numbers}
1
2
3
4
IBRION=-1
NSW=0
LWAVE = .T. (Write WAVECAR )
LCHARG = .T. (Write CHGCAR )

这一步完成后的输出文件CHGCAR即是电荷密度,导入VESTA直接查看。

bader电荷

自洽时INCAR再加一个参数:

{.line-numbers}
1
2
3
IBRION=-1
NSW=0
LAECHG=.T.

输出文件:
AECCAR0 全电子中的芯电子部分(POTCAR)
AECCAR1 自洽计算前的交叠原子电荷(POTCAR,bader 分析不用)
AECCAR2 全电子中的价电子部分(自洽计算后)
CHGCAR 自洽计算后赝电子

输入:

{.line-numbers}
1
2
chgsum.pl AECCAR0 AECCAR2
bader CHGCAR -ref CHGCAR_sum

得到两个文件:

  • BCF.dat 电荷极大值序号和坐标
    体积内的 Bader 电荷(积分)
    距离极大值最近的原子和距离
  • ACF.dat 原子序号和坐标
    原子电荷(极大值的电荷之后)
    到零通量面最小距离
    原子体积

李强博士《learn vasp the hard way》https://www.bigbrosci.com/2011/12/24/A08/ Bader电荷与原子电荷着色:
https://www.bigbrosci.com/2011/12/23/A07/ https://mp.weixin.qq.com/s/BNzjhz8SI_HXkaEnyL29Bg

成键反键

自洽时添加:

{.line-numbers}
1
2
3
4
5
NBANDS=96
ISYM=-1
IBRION=-1
NSW=0
LWAVE = T

vi lobsterin
原子编号改为自己的,例如

lobsterin:

1
2
3
4
5
6
7
8
COHPstartEnergy  -14
COHPendEnergy 6
basisSet Bunge
basisfunctions Ga 4s 4p 3d
basisfunctions As 4s 4p

cohpbetween atom 1 and atom 2

前两行代表提取的能量范围;第三行代表选择的基组,如果不设置(空行)就是默认基组,如果要覆盖默认基集,可使用此关键字进行覆盖。它目前支持bunge、koga和pbeVaspFit2015;第四行代表为每个元素指定所使用的基函数;最后一行代表选择一对两个原子进行键合分析。

lobster < lobsterin

1、删除存在的WAVECAR文件,做单点计算,计算需要保存波函数文件;
2、做完单点计算在当前目录下运行lobster软件;
3、用wxDragon软件打开COHPCAR.lobster文件。
4、COHPCAR.lobster 文件包含COHP结果,第二列为pCOHP,将第二列乘-1可得-pCOHP
https://github.com/JaGeo/LobsterPy
https://zhuanlan.zhihu.com/p/470592188

注:

  • orbitalWise 这个参数一般加在cohpBetween或者cohpGenerator之后,用以把原子之间的相互作用具体分到s p d 这些轨道上,比如计算s-s, p-p, d-d之间的pCOOP, pCOHP;
    比如:cohpBetween atom 1 atom 2 orbitalWise

cohpGenerator from 1.4 to 1.5 type Ga type Sr orbitalWise

  • 不要使用US-PP超软赝势,使用PAW赝势;请使用vasp_std版本,不要使用vasp_gam
    ;不支持处理对称:ISYM=0或ISYM=-1

  • 由于需要把基于平面波的波函数投影到原子轨道的基函数上,所以计算的能带的数量要多于投影的基函数,VASP计算的时候,默认的NBANDS达不到要求,需要加大NBANDS,只要查看OUTCAR里面的NBANDS多少,grep NBANDS OUTCAR;在INCAR里面,NBANDS关键词设置成1.5倍NBANDS,足以满足计算的投影要求,在计算过程中一定要记住设置足够的NBANDS,否则会出现报错,一直不能正常算出COHP

  • cohpGenerator 这里通过距离或者元素类别告诉想要计算什么原子之间的COHP;cohpGenerator from 1.5 to 2.3 type Fe type N orbitalwise #输出Fe-N原子对距离在1.5到2.3的所有数据之和。

  • gaussianSmearingWidth 高斯展宽,如果VASP中用的ISMEAR=0,那么gaussianSmearingWidth的数值和SIGMA一样,默认值是0.2 eV,但是过大了,所以需要减小。设置格式:gaussianSmearingWidth 0.05

ELF电子局域函数计算

电子局域函数(Electron Localization Function—ELF)是研究电子结构的手段之一。用来描述以某个位置处的电子为参考,在其附近找到与他同自旋的电子的概率,可以表征这个作为参考的电子的局域化程度,也是一种描述在多电子体系中的电子对概率的方法。
ELF的取值范围是[0,1]。
ELF=1 对应完全局域化;
ELF=1/2,对应类电子气型的成对概率;
ELF=0,对应电子完全离域化。

计算流程:做自洽的时候在INCAR中加入:

{.line-numbers}
1
2
3
4
PREC=A
LELF=T
IBRION=-1
NSW=0

运行vasp得到ELFCAR文件。

将ELFCAR拖入Vesta http://www.jp-minerals.org/vesta/en/download.html

点击 Utiltties - 2D Data Display。
点击Slice,确定切面和fractional coordinates。
导出图片:之后可以用PS工具将原子模型加入。

费米面

费米面计算,这里我们以Cu为例,见vaspkit/examples/Cu_fermi_surface。
第一步,准备好POSCAR,注意一定是原胞(primitive cell )和用于静态计算的INCAR文件;
第二步,运行VASPKIT,输入261命令,产生用于计算计算费米面的KPOINTS和POTCAR文件,
第三步,提交VASP作业;第四步,待VASP计算完成后,再次运行VASPKIT,输入262命令,生成包含费米面数据的FERMISURFACE.bxsf文件;第五步,运行XcrySDen软件。
发现第五条能带穿过费米面,因此我们选择Band Number 5,然后点击Selected,然后得到Cu的费米面。推荐使用FermiSurfer可视化费米面。
对于半导体不存在费米面,如果想可视化导带或价带中某个能级处的等能面,新建FERMI_ENERGY.in文件然后在该文件第二行(第一行为注释行)中人为指定能级大小,注意该值一定要保证穿过价带或导带。接下来调用vaspkit-252或263命令提取可视化。

STM模拟

参考:https://zhuanlan.zhihu.com/p/538950799

https://zhuanlan.zhihu.com/p/667558468

https://www.vasp.at/wiki/index.php/STM_of_graphite

  1. 首先做一个自洽计算,输出 WAVECAR 和 CHGCAR 文件。

  2. 加入一些特殊参数计算,输出我们想要的 PARCHG 文件。

    1
    2
    3
    4
    5
    LPARD = .TRUE.
    LSEPK = .FALSE.
    LSEPB = .FALSE.
    NBMOD = -3
    EINT = -0.1 0.1
  3. 用p4vasp软件打开PARCHG文件或者vaspkit 后处理。由于 VASP 在对于远离原子表面时波函数的衰减指数描述不明确,所以模拟探针高度的时候不应该太高,否则会出现非物理的现象,比如 7 埃,就算比较远了。所以一般高度的选择要接近材料表面,这里选择了 0.5 埃。而原胞沿着 ab 平面重复的标准是 20 埃,也就是说原胞晶格 * 重复单元 = 20 埃就可以。vaspkit 里选择的是 Tersoff-Hamman 方法,现在大多数模拟 STM 的都是选择这个方法。

  4. vaspkit处理后的输出文件包括X.grd, Y.grd, STM_*.grd三个文件,可以用以下脚本画图。
    Usage: python plot.pyimport numpy as npimport matplotlib as mpl#mpl.use(‘Agg’) #silent modefrom matplotlib import pyplot as plt
    x_mesh=np.loadtxt(‘X.grd’) y_mesh=np.loadtxt(‘Y.grd’) v_mesh=np.loadtxt(‘STM.grd’)
    fig = plt.figure()cmap = plt.cm.get_cmap(“Greys_r”) plt.contourf(x_mesh,y_mesh,v_mesh,cmap=cmap)
    plt.axis(‘equal’)plt.axis(‘off’)#plt.show()plt.savefig(“STM.jpg”,dpi=300)
    或者用以下脚本处理后,得到data.txt文件,将该文件导入origin画图。Usage: sh get_data.shawk ‘{for(i=1;i<=NF;i++) print $i}’ X.grd >xawk ‘{for(i=1;i<=NF;i++) print $i}’ Y.grd >yawk ‘{for(i=1;i<=NF;i++) print $i}’ STM_0.50.grd >spaste x y s >data.txt

态密度

自洽计算之后,读取产生的CHGCAR、WAVECAR
上面例如bader电荷、ELF标签都关闭。
INCAR中修改:

{.line-numbers}
1
2
3
4
5
ISTART=1
ICHARG=11
ISMEAR=-5
NEDOS=2001
LORBIT=11

其他参数和自洽时保持一致。
运行后vaspkit的111/112/113/114命令处理即可得到。

分波态密度和局域态密度

态密度很多分析和能带的分析结果可以一一对应,很多术语也和能带分析相通。但是因为它更直观,因此在结果讨论中用得比能带分析更广泛一些。在电子能级为准连续分布的情况下,单位能量间隔内的电子态数目,即为态密度。能态密度与能带结构密切相关,是一个重要的基本函数。固体的许多特性,如电子比热、光和X射线的吸收和发射等,都与能态密度有关。在VASP中,LORBIT=10计算的就是LDOS,也就是每个原子的局域态密度 (local DOS),是分解到每个原子上面的s,p,d;LORBIT=11,计算的就是PDOS,投影态密度(projected DOS)或分波态密度(partial DOS),不仅分解到每个原子的s,p,d,而且还进行px,py,pz分解。

/vaspkit.0.72/examples/LDOS_PDOS/Partial_DOS_of_CO_on_Ni_111_surface 目录下启动 vaspkit ,输入命令 115 选择子功能 The Sum of Projected DOS for Selected Atoms and orbitals 。我们需要提取的是O的s,p轨道,C的s,p轨道以及表面的dxy和d3z2-r2(dz2)轨道。首先提示你选择元素(累加),第一次输入元素O,回车后提示你输入提取的轨道名(累加),第一次输入轨道s。回车后重复以上两个操作。如果想结束输入,在元素选择行直接按回车键结束输入。

元素行接受自由格式输入,你可以输入以下内容元素行接受自由格式输入,你可以输入以下内容1-3 4 Ni表示选择元素1,2,3,4和Ni元素的PDOS进行累加。轨道行只支持标准输入,如果使用了LORBIT=10不投影轨道,你只能选择 s p d f中的一个或多个,如果使用了LORBIT=11投影了轨道,你可以从s p d f和s py pz px dxy dyz dz2 dxz dx2 f-3 f-2 f-1 f0 f1 f2 f3中选择一个或多个输入。轨道行还支持输入all计算所有轨道的DOS之和。

比如元素行输入C O ,对应的轨道选择s px,那么提取的就是所有C和O元素的s轨道和px轨道的PDOS之和。 PDOS_USER.dat中的轨道名 为#Energy C&O_s&px,显而易懂。

电荷密度差

Vaspkit 可以方便的处理 VASP 输出的电荷密度文件 (CHGCAR),做电荷密度差分处理,VTST 脚本或者 VESTA 也可以做,但是不如 vaspkit方便灵活。
注:此文中讲述的功能同样适用于任何的实空间函数形式的文件处理,比如静电势差值图 (LOCPOT),电子定域化函数差值图 (ELFCAR),能带分解的电荷密度插值图 (PARCHG) 等等。
电荷密度差分(charge density difference)是研究电子结构的重要手段之一。可以直观的得到两个片段相互作用后的电子流向,原子形成分子过程中电子密度的变化、探究化学键的本质。电荷密度差分有以下几种形式:
(1)体系的电荷密度减去组成它的两个或几个片段的密度,也是最常见的电荷密度差分形式:

(2)自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度),也称为变形电荷密度(Deformation charge density)

(3)在某个状态的密度减去这个体系在另外一个状态的密度。比如:外加电场作用下的电荷密度减去没有外势场的电荷密度。再比如:激发态的密度减去它在基态时的密度。

以上三种电荷密度差,无论哪一种计算都需要保证原子坐标坐标必须一致!并且保证晶胞参数和 FFT-mesh 的格点密度 (NGXF, NGYF, NGZF) 一致!

(1) 的计算只需要优化 AB 的结构,A、B片段保持 AB 结构中原子坐标,不能再对片段进行单独优化。
(2) 只要优化自洽计算时候的几何结构,生成原子的球对称的电荷密度时候 (ICHARG = 12) 保持坐标不变。
(3) 也需要共用其中一个状态的几何构型(比如激发态或基态的几何构型)。

CHGCAR文件格式

CHGCAR 是包含电子密度信息的格点文件,对于自旋非极化体系(ISPIN = 1)计算只包含电荷密度,对于自旋极化体系(ISPIN = 2)计算还包含自旋电子密度。可以用 VESTA,Jmol 等程序打开。CHGCAR 的第一部分和 POSCAR, CONTCAR 的格式是完全一样的,它包含了最终结构的晶格矢量,原子核坐标等信息。紧接着是三个数字,这三个数是实空间函数的网格密度,对应于 NGXF,NGYF,NGZF 三个变量。然后是电荷密度信息 ρ(r) * Vcell。一共有 NGXF * NGYF * NGZF 个数值。比如:下面是一个 Al2O3 晶胞的 CHGCAR 文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Al2O3 Cell opt                          
1.00000000000000
6.138452 0.000551 -0.011568
5.434311 2.854620 -0.011568
-1.324042 -0.326661 5.501689
O Al
6 4
Direct
0.160314 0.160314 0.108934
0.839686 0.839686 0.891066
0.495028 0.495028 0.257399
0.504972 0.504972 0.742601
0.826429 0.826429 0.432705
0.173571 0.173571 0.567295
0.090130 0.090130 0.795615
0.909870 0.909870 0.204385
0.342109 0.342109 0.683236
0.657891 0.657891 0.316764
100 100 96
0.43432567054E+01 0.44027037758E+01 0.45822813500E+01 0.48861003640E+01 0.53223812855E+01
0.59054783755E+01 0.66583226816E+01 0.76145085758E+01 0.88190771181E+01 0.10327304760E+02
...

Vaspkit计算几个片段的电荷密度差

以 CO 分子吸附在 Ni(100) 表面为例:
步骤一:优化 CO/Ni(100) 的结构;
步骤二:分别计算 Ni(100) 和 CO 的单点能 CO 和 Ni(100) 片段的坐标从 CO/Ni(100) 的 CONTCAR 里直接截取,不要再结构优化!!计算时也要保证三次自洽计算所采用的 FFT mesh 一致(NGXF,NGYF,NGZF)
步骤三:用 vaspkit 314 功能做电荷差分。依次输入三个片段的 CHGCAR 路径:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
314   
======== File Options ========
Input the Names of Charge/Potential Files with Space:
(e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR)
(e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)

------------>>
./CHGCAR ./co/CHGCAR ./slab/CHGCAR

-->> (01) Reading Structural Parameters from ./CHGCAR File...
-->> (02) Reading Charge Density From ./CHGCAR File...
-->> (03) Reading Structural Parameters from ./co/CHGCAR File...
-->> (04) Reading Charge Density From ./co/CHGCAR File...
-->> (05) Reading Structural Parameters from ./slab/CHGCAR File...
-->> (06) Reading Charge Density From ./slab/CHGCAR File...
-->> (07) Written CHGDIFF.vasp File!

生成 CHGDIFF.vasp 包含电荷密度差的信息,可以直接导入到 VESTA 里作图。
图片

Vaspkit计算变形电荷密度差

变形电荷密度差:自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度)。
以 CO 分子为例:
步骤一:先自洽计算优化 CO 分子。
步骤二:新键文件夹,在优化好的结构基础上用 ICHARG = 12 做非自洽计算。
步骤三:用 vaspkit 314 功能做电荷差分。依次输入两个 CHGCAR 路径。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
314   
======================= File Options ============================
Input the Names of Charge/Potential Files with Space:
(e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR)
(e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)

------------>>
./CHGCAR ./deform/CHGCAR

-->> (01) Reading Structural Parameters from ./CHGCAR File...
-->> (02) Reading Charge Density From ./CHGCAR File...
-->> (03) Reading Structural Parameters from ./deform/CHGCAR File...
-->> (04) Reading Charge Density From ./deform/CHGCAR File...
-->> (05) Written CHGDIFF.vasp File!

步骤四:得到 CHGDIFF.vasp 文件,导入到 VESTA 里做图。

Vaspkit计算外加电场下的电荷密度差

外加电场作用下的电荷密度减去没有外势场的电荷密度。
以 InSe 二维单层材料为例:
步骤一:先优化没有外加电场的结构。
步骤二:在同样结构下计算外加电场下做单点自洽计算。添加关键词,EFIELD 控制电场力的大小(eV/Angstrom)F=qE, 对应电场的单位是 E=F/q, 场强单位是 V/Angstrom

1
2
3
EFIELD = 0.05
IDIPOL = 3
LDIPOL = .TRUE.

步骤三:运行 Vaspkit 314 功能做电荷差分。依次输入二个 CHGCAR 路径。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
314
======================= File Options ============================
Input the Names of Charge/Potential Files with Space:
(e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR)
(e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)

------------>>
./CHGCAR ./noele/CHGCAR

-->> (01) Reading Structural Parameters from ./CHGCAR File...
-->> (02) Reading Charge Density From ./CHGCAR File...
-->> (03) Reading Structural Parameters from ./noele/CHGCAR File...
-->> (04) Reading Charge Density From ./noele/CHGCAR File...
-->> (05) Written CHGDIFF.vasp File!
步骤四:得到 CHGDIFF.vasp 文件,导入到 VESTA 里做图。

电荷密度差作图

得到的电荷密度差文件必须作成图才便于被直观地考察,常用的图分为三类:
(1)3D 等值面图,isosurface.
(2)2D 切面图.
(3)1D,平面平均的数值,planar averaged charge.
通常这几种图可以配合使用讨论。等值面图和切面图用 VESTA 画最方便,平面平均图用 vaspkit 生成方便。得到的 CHGDIFF.vasp 用 VESTA 打开,isosurface 青色部分电荷密度减小,黄色密度电荷密度增加。注意 VESTA 对电荷密度的默认单位是 e/Bohr3。
vaspkit 319 功能可以把VASP的格点文件转化成 .cube的方法,.cube 在 VMD 里作图可以得到更好看的效果:

能带计算

参考链接:
https://yh-phys.github.io/2019/10/04/vasp-band-structure-dos/
自洽计算之后,读取产生的CHGCAR、WAVECAR

{.line-numbers}
1
2
ISTART=1
ICHARG=11

其他参数和自洽时保持一致。

能带计算的KPOINTS文件不同于前面,称为高对称点。
学术之友微信公众号历史记录搜索文章《五种方法生成能带结构计算的高对称点》
1、http://cryst.ehu.es/
2、https://www.materialscloud.org/work/tools/seekpath
3、http://www.aflowlib.org/aflow-online (已更新)
4、http://pymatgen.org/ (pymatgen)
5、http://www.xcrysden.org/ (XCrySDen)
Accurate DOS and Band-structure calculations:
http://cms.mpi.univie.ac.at/vasp/vasp/Accurate_DOS_Band_structure_calculations.html

可以通过vaspkit的303命令产生KPATH.in文件,cp KPATH.in KPOINTS即可。
计算完成后通过vaspkit的211命令提取能带。
注意能带路径生成操作需在自洽计算前完成。

  • 1.准备POSCAR,调用vaspkit- 303(体相材料)或-302(二维材料)得到KPATH.in和PRIMCELL.in文件;对于二维体系,需要检查POSCAR文件的真空层是否沿z方向,如果没有,可调用vaspkit-923或vaspkit-407强制真空层沿z方向。另外,如果需要改变识别结构对称性精度(symprec默认值为1E-5),可通过vaspkit -task 302 -symprec 1E-6实现或者在~/.vaspkit中修改SYMPREC默认值。

  • 2.执行cp PRIMCELL.vasp POSCAR命名。(如果PRIMCELL.vasp 和POSCAR不一样,并且要读取电荷密度和波函数,要从SCF计算保持结构一致。)KPATH.in推荐能带路径只针对于PRIMITIVE CELL,缺少这一步,你可能得到错误的结果。如果有必要,比较KPATH.in文件中的能带路径是否与在线能带路径产生工具SeeK-Path产生的一致,包括比较PRIMCELL.vasp和HIGH_SYMMETRY_POINTS文件。需要指出的是SeeK-Path只用于体相结构能带路径的产生。自1.3.0版本,对于二维体系不需要执行cp PRIMCELL.vasp POSCAR。

投影电荷密度

复制能带结构文件,
vi INCAR

{.line-numbers}
1
2
3
增加标签 LPARD=.T.
增加标签 IBAND=16 17
增加标签 LSEPB=.T.

  • LPARD:它的取值为.TRUE.或.FALSE.,它的默认值是.FALSE.,当为.TRUE.时,表示读入自
    洽收敛的 CHGCAR 和 WAVECAR 并进行 partial charge density计算。
  • IBAND:它的取值就是你想要计算的第几条能带或哪几条能带(比如要计算第 4、 5、 6 条能带,那么就设置 IBAND = 4 5 6),此时 NBMOD的值就是所要计算的能带的条数。
  • EINT:另一种方式指定所要计算的能带,它是指定计算某能量范围的partial charge density,一般是设置为两个实数,比如EINT= 4.00 5.00,此时NBOMD的值设置为-2。如果只设置了一个数,那么表示计算从 EINT 到费米能级这个范围内的,此时NBOMD 的值为-3。
  • KPUSE:指定所要计算的 k 点(哪个或哪几个)。
  • LSEPB:指定是不是要把计算的partial charge density按每个带分别写到各自对应的文件
  • PARCHG.nb (设置为.TRUE.)中,还是把它们合并写到一个文件中(相当于把各个带对应的partial charge density 加起来,设置为.FALSE.),默认值为.FALSE.。
  • LSEPK:指定是不是把要计算的partial charge density按每个k分别写到各自对应的文件PARCHG.nk (设置为.TRUE.).

三个方法:

1、读取自恰计算的波函数,KPOINTS文件采用特定的k点

{.line-numbers}
1
2
3
4
5
6
#partial charge density
LPARD = .TRUE.
IBAND = 17 18 # 价带顶和导带底
KPUSE = 1
LSEPB = .TRUE.
LSEPK = .TRUE.
{.line-numbers}
1
2
3
4
CsPbI3-cubic
1
r
0.500000 0.500000 0.500000 1.0

http://muchong.com/html/200906/1404880.html
2、读取自恰计算的波函数,KPOINTS文件采用自恰计算的k点

{.line-numbers}
1
2
3
4
5
6
7
8
#partial charge density
LPARD = .TRUE.
IBAND = 17 18
KPUSE = 35
LSEPB = .TRUE.
LSEPK = .TRUE.
NELM = 200
ICHARG = 1

在scf自恰计算OUTCAR中
找到对应导带底和价带顶
对应的k-point编号
{.line-numbers}
1
2
3
4
5
CsPbI3-cubic
0
Gamma
8 8 8
0.0 0.0 0.0

3、读取能带计算的波函数,KPOINTS文件采用能带计算的k点

{.line-numbers}
1
2
3
4
5
6
7
8
#partial charge density
LPARD = .TRUE.
IBAND = 17 18
KPUSE = 80 81
LSEPB = .TRUE.
LSEPK = .TRUE.
NELM = 200
ICHARG = 1

KPOINTS和能带计算一样。在band能带计算OUTCAR中找到对应导带底和价带顶对应的k-point编号。

轨道投影态密度和能带。

{.line-numbers}
1
2
LORBIT=11 
NEDOS=2001

提取DOS:vaspkit -task 111/113
提取band:211/212/213/214
如果3D能带图不够平滑,可通过能带插值进行改善。
VASPKIT升级到1.3.0或更新版本,在~/.vaspkit文件中设置以下三个参数:

{.line-numbers}
1
2
3
GET_INTERPOLATED_DATA      .TRUE.     是否插值
INTERPOLATION_SPACING 0.04 插值间距,决定平滑程度
INTERPOLATION_METHOD 'cubic' 插值方法

HSE06计算电子结构:

  1. 结构优化(PBE)

  2. 静态自洽scf计算(PBE)KPOINTS和POTCAR文件不变,写波函数和电荷密度

  3. (可选)静态自洽scf计算(HSE06)KPOINTS,POSCAR和POTCAR文件不变

  4. DOS/band计算(HSE06)vaspkit -task 303
    vaspkit→251→2→0.03→0.03生成KPOINTS,运行VASP
    vaspkit-252
    http://vaspkit.cn/index.php/29.html

  5. PBE优化结构

  6. HSE自洽计算(得到波函数)

  7. HSE能带计算/态密度计算

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    ALGO = Damped
    HFSCREEN = 0.2
    LHFCALC = .True.
    LMAXFOCK = 4
    PRECFOCK = Fast
    TIME = 0.4

+SOC计算电子结构:

  1. 结构优化(PBE)
  2. 静态自洽scf计算(PBE+SOC)KPOINTS和POTCAR文件不变,写波函数和电荷密度
  3. 非自洽计算(PBE+SOC)
    Spin-Orbit Coupling Calculation
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
     ISPIN  =  2
    NELMIN = 6 (Min electronic SCF steps)
    LSORBIT = .TRUE. (Activate SOC)
    GGA_COMPAT = .FALSE. (Apply spherical cutoff on gradient field)
    #VOSKOWN = 1 (Enhances the magnetic moments and the magnetic energies)
    LMAXMIX = 2 (For d elements increase LMAXMIX to 4, f: LMAXMIX = 6)
    SAXIS = 0 0 1 (Direction of the magnetic field)
    #MAGMOM = 0 0 3 (Set this parameters manually, Local magnetic moment parallel to SAXIS, 3*NIONS*1.0 for non-collinear magnetic systems)
    NBANDS = 32 (Set this parameters manually, 2 * number of bands of collinear-run)

调节SOC强度

  1. 调节自旋轨道耦合作用的大小

要调优VASP中的自旋轨道耦合强度,需要修改源文件vasp_source_code_path/src/relativistic.F 并重新编译。如果要将强度减少一半,可以在relativistic.F的129行L●S项乘以0.5d0,如下:
DO I=0,1DO J=0,1DO M =1,2LL+1DO MP=1,2LL+1    DLLMM(LMP+MP-1,LM+M-1,J+2I+1)=DLLMM(LMP+MP-1,LM+M-1,J+2I+1)+ &    0.5d0SUMLS(M,MP,I+2*J+1,LL)   !!!line 129 relativistic.F fileEND DOEND DOEND DOEND DO

代码片段:可切换语言,无法单独设置文字格式

   注意,需要使用重新编译的vasp_ncl来运行计算。另外,最好比较0d0L●S,1d0L●S和vasp_std,vasp_ncl计算结果。

  1. 如果只想给一些元素加SOC,而其他元素不加该怎么办呢?拓扑材料计算互助群的Eurus-清华给出了VASP中给特定元素加自旋轨道耦合相互作用的方法(感谢Eurus-清华)。在原来VASP,src文件夹下的relativistic.F文件中注释掉46行的

! REAL(q), PARAMETER :: INVMC2=7.45596E-6

然后53行位置加入如下6-17行的内容,然后重新编译。

  INTEGER, PARAMETER :: LMAX=3, MMAX=LMAX*2+1      COMPLEX(q) DUMMY(MMAX,MMAX,3,LMAX)      COMPLEX(q) LS(MMAX,MMAX,4,LMAX)      REAL(q) SUM, SCALE

! add  parameter  REAL(q), PARAMETER :: INVMC2_ORIG=7.45596E-6 REAL(q) INVMC2! add end ! add  control IF(ABS(Z-83).LT.1.0E-5)THEN INVMC2=INVMC2_ORIG ELSE INVMC2=INVMC2_ORIG*0.0001 END IF! add end
LS=(0._q,0._q) CALL SETUP_LS(1,THETA,PHI,DUMMY(1:3,1:3,1:3,1),LS(1:3,1:3,1:4,1)) CALL SETUP_LS(2,THETA,PHI,DUMMY(1:5,1:5,1:3,2),LS(1:5,1:5,1:4,2)) CALL SETUP_LS(3,THETA,PHI,DUMMY(1:7,1:7,1:3,3),LS(1:7,1:7,1:4,3))
原理:这样做的原理就是将除Z-83的元素外的其他元素缩小一万倍。注意:这个脚本的第12行, IF(ABS(Z-83)),表明除了83号元素打开SOC之外,其他所有的元素SOC都是相当于关闭的。如果想打开某个特定元素的SOC,就把83改成你想要的元素序号就行。
另外,其他元素SOC关闭方式就是将其缩小一万倍。所以使用这种方法,任然是非共线计算,如果要进行wannier90计算,只要按照正常的wannier90设置就可以。

HSE06+SOC band structure calculation

  • Step-I: HSE06 SCF calculation with SOC (INCAR):

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    #SOC-related:
    LSORBIT=.TRUE.
    SAXIS= 0 0 1
    ISYM=0

    #HSE06-related:
    LHFCALC = .TRUE.
    HFSCREEN = 0.2
    ALGO = Damped #Damped/Normal/All
    TIME = 0.4
  • Step-II: HSE06 NSCF calculation with SOC:

    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    #SOC-related:
    LSORBIT=.TRUE.
    SAXIS= 0 0 1
    ISYM=0

    #HSE06 related:
    LHFCALC = .TRUE.
    HFSCREEN = 0.2
    ALGO = Normal !use Normal
    TIME = 0.4

能带反折叠

vaspkit:

  1. 准备POSCAR
  2. 运行vaspkit -task 302/303生成KPATH.in文件
  3. cp PRIMCELL.vasp POSCAR
  4. 运行vasp优化结构,cp CONTCAR POSCAR
  5. 运行vaspkit -task 400/400生成超胞结构POSCAR,并执行cp TRANSMAT TRANSMAT.in
  6. 运行vasp优化含有缺陷的超胞结构,cp CONTCAR POSCAR
  7. 运行vaspkit -task 281生成KPOINTS文件
  8. 准备INCAR(静态计算,设置LWAVE = T,适当增加NBANDS值)并执行VASP计算
  9. 运行vaspkit -task 282提取effective band structure(画法类似于投影能带)
  • 如果不会计算变换矩阵(TRANSMAT.in文件不存在时),可以执行cp PRIMCELL.vasp PRIMCELL.in,281 ,这时vaspkit分别读入超胞基矢(POSCAR)和原胞基矢(PRIMCELL.in)计算变换矩阵并生成TRANSMAT.in文件。
  • vaspkit/examples/band_unfolding/ebs.py画图。生成的ENERGY.grd,WEIGHT.grd和MOMENTUM.grd三个文件,可使用vaspkit/examples/band_unfolding/ebs_k_resolved_dos_plot_matlalb.m脚本画图。

bandup:

第一步 得到超胞CHGCAR:输入文件参考BANDUP的example的第一步,提交任务之后得到CHGCAR。
第二步 主要是为了得到超胞非自洽计算的K点路径,分别给出单胞和超胞的晶格坐标文件prim_cell_lattice.in与supercell_lattice.in。还有单胞的高对称K点路径。
下一步命令是bandup kpts-sc-get -no_symm_avg,这一步会生成KPOINTS_supercell.out文件后面会用到。-no_symm_avg这个标签是为了减少需要计算的k点个数,减小计算量和WAVECAR大小。
第三步 利用前面第一步生成 CHGCAR; POSCAR POTCAR与第一步的一样,K点从第二步那个复制过来,然后提交任务得到波函数WAVECAR
第四步就要开始做能带反折叠了,一句命令,bandup unfold -emin -13 -emax 6 -dE 0.050 -no_symm_avg
会生成一个unfolded_EBS_not-symmetry_averaged.dat文件,画出来就行了。
https://github.com/band-unfolding/bandup

http://staff.ustc.edu.cn/~zqj/posts/Band-unfolding-tutorial/
可以处理能带反折叠的软件汇总:
https://mp.weixin.qq.com/s/-5lN2lwxsALA-rhNQxK1dw
https://vaspkit.com/tutorials.html#band-unfolding
https://github.com/QuantumMaterialsModelling/bands4vasp
https://github.com/QuantumMaterialsModelling/UnfoldingPatch4vasp

3D能带

__第一步__,准备好石墨烯的POSCAR和用于静态计算的INCAR文件;
__第二步__:运行VASPKIT并选择231生成用于3D能带计算 的KPOINTS文件。为了获得平滑的3D能带面,用于产生KPOINTS文件的分辨率值应该设置在0.005左右 ;
__第三步__,提交VASP作业;
__第四步__,待VASP计算完成后,再次运行VASPKIT,输入232或233命令。233命令可一次性输入包含价带顶的BAND.HOMO.grd和导带底的BAND.LUMO.grd文件的能带;232可以得到其它任意能带的计算数据,这里我们选择233。
__第五步__:运行python how_to_visual.py,绘制3D能带图(python2.7环境)。
如果3D能带图不够平滑,可通过能带插值进行改善。

有效质量:

除了通过能带,也可以通过调用VASPKIT-911命令来得到带边位置,接下来我们计算K点处沿着K -> Γ和K -> M方向的电子及空穴载流子的有效质量;
第一步:准备POSCAR文件以及VPKIT.in文件,VPKIT.in文件包含以下内容;

{.line-numbers}
1
2
3
4
5
6
1                  设置为1将会生成计算有效质量的KPOINTS文件,2则计算有效质量
6 6表示我们在K附近取6个离散点用于多项式拟合
0.015 k-cutoff, 截断半径,典型值0.015 1
2 2 表示有两个有效质量任务
0.333333 0.3333333 0.000 0.000 0.000 0.000 K->Γ 计算K点处有效质量,沿着K指向Γ方向
0.333333 0.3333333 0.000 0.500 0.000 0.000 K->M 计算K点处有效质量,沿着K指向M方向

第二步:运行vaspkit并选择912或者913产生KPOINTS,POTCAR文件及静态计算INCAR文件;
第三步:提交VASP作业;
第四步:计算完成后把VPKIT.in文件中第一行的1修改为2,然后再次运行vaspkit并选择913,得到结果:
为了得到可靠的有效质量计算值,括号里的精度至少要在1E-7这个量级,如果误差比较大,可以适当调整k-cutoff值,一般减少该值但不能太小。
手动vaspkit处理的能带:$m^∗=3.815/B$

GW方法

参考:https://www.vasp.at/wiki/index.php/Bandgap_of_Si_in_GW

步骤:

  1. PBE基态计算
  2. 获得虚拟轨道(50-100/atom)

INCAR:

1
2
3
4
5
6
7
8
9
10
11
ALGO = Exact
NBANDS = 64
LOPTICS = .TRUE. #得到WAVEDER文件
CSHIFT = 0.1
NEDOS = 2000
# you might try
#LPEAD = .TRUE.

ISMEAR = 0
SIGMA = 0.05
EDIFF = 1E-8

完成后:

1
2
cp WAVECAR WAVECAR.DIAG
cp WAVEDER WAVEDER.DIAG
  1. GW计算
    读取上一步
    INCAR:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    # Frequency dependent dielectric tensor including
    # local field effects within the RPA (default) or
    # including changes in the DFT xc-potential (LRPA=.FALSE.).
    # N.B.: beware one first has to have done a
    # calculation with ALGO=Exact, LOPTICS=.TRUE.
    # and a reasonable number of virtual states (see above)
    ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50

    # be sure to take the same number of bands as for
    # the LOPTICS=.TRUE. calculation, otherwise the
    # WAVEDER file is not read correctly
    NBANDS = 64

    # Add this to update the quasiparticle energies
    # in the Green's function (GW0)
    #NELM = 4

    ISMEAR = 0
    SIGMA = 0.05
    EDIFF = 1E-8
    在OUTCAR中quasi-particle (QP) energies.
    1
    2
    3
    4
    5
    6
    7
    QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
    for sc-GW calculations column KS-energies equals QP-energies in previous step
    and V_xc(KS)= KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1 + <V_x>^1)

    k-point 1 : 0.0000 0.0000 0.0000
    band No. KS-energies QP-energies sigma(KS) V_xc(KS) V^pw_x(r,r') Z occupation
    .........

LWANNIER90计算能带结构

重新开始第三步

INCAR:

1
2
3
4
5
6
7
8
ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50
#LRPA = .FALSE.
## be sure to take the same number of bands as for
## the LOPTICS=.TRUE. calculation, otherwise the
## WAVEDER file is not read correctly
NBANDS = 64
##VASP2WANNIER90
LWANNIER90=.TRUE.

wannier90.win:

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
num_wann=8
num_bands=8

exclude_bands = 9-64

Begin Projections
Si:sp3
End Projections

dis_froz_max=9
dis_num_iter=1000

guiding_centres=true

# Bandstructure plot
#restart = plot
#bands_plot = true
#begin kpoint_path
#L 0.50000 0.50000 0.5000 G 0.00000 0.00000 0.0000
#G 0.00000 0.00000 0.0000 X 0.50000 0.00000 0.5000
#X 0.50000 0.00000 0.5000 K 0.37500 -0.37500 0.0000
#K 0.37500 -0.37500 0.0000 G 0.00000 0.00000 0.0000
#end kpoint_path
#bands_num_points 40
#bands_plot_format gnuplot xmgrace

begin unit_cell_cart
2.7150000 2.7150000 0.0000000
0.0000000 2.7150000 2.7150000
2.7150000 0.0000000 2.7150000
end unit_cell_cart

begin atoms_cart
Si 0.0000000 0.0000000 0.0000000
Si 1.3575000 1.3575000 1.3575000
end atoms_cart

mp_grid = 4 4 4

begin kpoints
0.0000000 0.0000000 0.0000000
0.2500000 0.0000000 0.0000000
0.5000000 0.0000000 0.0000000
0.2500000 0.2500000 0.0000000
0.5000000 0.2500000 0.0000000
-0.2500000 0.2500000 0.0000000
0.5000000 0.5000000 0.0000000
-0.2500000 0.5000000 0.2500000
0.0000000 0.2500000 0.0000000
0.0000000 0.0000000 0.2500000
-0.2500000 -0.2500000 -0.2500000
-0.2500000 0.0000000 0.0000000
0.0000000 -0.2500000 0.0000000
0.0000000 0.0000000 -0.2500000
0.2500000 0.2500000 0.2500000
0.0000000 0.5000000 0.0000000
0.0000000 0.0000000 0.5000000
-0.5000000 -0.5000000 -0.5000000
0.0000000 0.2500000 0.2500000
0.2500000 0.0000000 0.2500000
-0.2500000 -0.2500000 0.0000000
-0.2500000 0.0000000 -0.2500000
0.0000000 -0.2500000 -0.2500000
0.0000000 0.5000000 0.2500000
0.2500000 0.0000000 0.5000000
-0.2500000 -0.2500000 0.2500000
-0.5000000 -0.2500000 -0.5000000
0.2500000 0.5000000 0.0000000
0.2500000 -0.2500000 -0.2500000
-0.5000000 -0.5000000 -0.2500000
0.0000000 0.2500000 0.5000000
-0.2500000 0.2500000 -0.2500000
-0.2500000 -0.5000000 -0.5000000
0.5000000 0.0000000 0.2500000
-0.5000000 -0.2500000 0.0000000
0.0000000 -0.5000000 -0.2500000
-0.2500000 0.0000000 -0.5000000
0.2500000 0.2500000 -0.2500000
0.5000000 0.2500000 0.5000000
-0.2500000 -0.5000000 0.0000000
-0.2500000 0.2500000 0.2500000
0.5000000 0.5000000 0.2500000
0.0000000 -0.2500000 -0.5000000
0.2500000 -0.2500000 0.2500000
0.2500000 0.5000000 0.5000000
-0.5000000 0.0000000 -0.2500000
0.0000000 -0.2500000 0.2500000
0.2500000 0.0000000 -0.2500000
-0.2500000 -0.2500000 -0.5000000
0.2500000 0.5000000 0.2500000
0.2500000 -0.2500000 0.0000000
-0.5000000 -0.2500000 -0.2500000
0.2500000 0.2500000 0.5000000
0.0000000 0.2500000 -0.2500000
-0.2500000 -0.5000000 -0.2500000
0.5000000 0.2500000 0.2500000
-0.2500000 0.0000000 0.2500000
0.0000000 0.5000000 0.5000000
0.5000000 0.0000000 0.5000000
0.2500000 -0.2500000 0.5000000
0.5000000 0.2500000 -0.2500000
-0.5000000 -0.2500000 -0.7500000
0.2500000 -0.5000000 -0.2500000
-0.2500000 0.2500000 -0.5000000
end kpoints

Compute Wannier functions
run wannier90:

wannier90.x wannier90

This run generates the wannier90 standard output (wannier90.wout) and the file wannier90.chk needed for the wannier interpolation (next step)

Obtain bandstructure (Wannier interpolation)
Uncomment the bandstructure plot flags in wannier90.win and rerun (restart) wannier90:

wannier90.x wannier90

This run generates the following bandstructure files which can be visualized using xmgrace or gnuplot:

wannier90_band.agr

wannier90_band.dat

wannier90_band.gnu

to plot the band structure using gnuplot: gnuplot -persist wannier90_band.gnu

Beyond the random-phase-approximation

To include local field effects beyond the random-phase-approximation in the description of the frequency dependent dielectric response function (local field effects in DFT) add the following line to your INCAR file:

LRPA = .FALSE.
and again restart from the WAVECAR and WAVEDER files from step 2.

Beyond G0W0: GW0

The most usual step beyond single-shot GW (G0W0) is to iterate the quasi-particle energies in the Greens functions. This is the socalled GW0 approximation. To have VASP do, for instance, 4 iterations of the QP-energies in G, add the following line to the INCAR file:

NELM = 4
and again restart from the WAVECAR and WAVEDER files from step 2.

To quickly find the QP-energy of the highest lying occupied state after 4 iterations of the QP energies in G, type:

./gap_GW.sh OUTCAR

截断能

截断能设置在INCAR文件的ENCUT参数。截断能(energy cutoff)指定了展开波函数用到的平面波基组(planewave basis set)所对应的截断能量。
参考: https://zhuanlan.zhihu.com/p/348826693
截断能越大,代表用来描述波函数的平面波基组越多,计算精度越高,同时计算耗时也越长。
如果没有指定ENCUT,程序将默认ENCUT=ENMAX(注意:当PREC = High时, ENCUT = 1.3*ENMAX)。扩胞并不会影响ENCUT的设置。

  • 测试标准
    ENMAX在保证定性正确的基础上大幅度降低计算耗时,如果想要定量计算某些性质就需要增大截断能了。所以ENCUT测试通常从ENMAX开始测试,一般截止到 1.5*ENMAX 即可(注意:比ENMAX小的测试可能会得到错误的能量)。当体系中每个原子的能量差收敛至0.001 eV/atom时,此ENCUT即为最终取值。

简单测试:
for i in {1..8}; do cp 400 $((400+$i*50)); sed -i "s/400/$((400+$i*50))/g" $((400+$i*50))/INCAR ; done

encut.sh脚本: 参考 https://zhuanlan.zhihu.com/p/115235135

{.line-numbers}
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#!/bin/bash

#生成KPOINTS文件,KPOINTS不能与INCAR中的KSPACING参数共用,二选一
cat > KPOINTS <<!
A
0
M
9 9 9
0 0 0
!

#产生计算所需POSCAR
cat > POSCAR<<!
Si8
1.0000000000
5.4687280000 0.0000000000 0.0000000000
0.0000000000 5.4687280000 0.0000000000
0.0000000000 0.0000000000 5.4687280000
Si
8
Direct
0.0000000000 0.0000000000 0.0000000000
0.2500000000 0.7500000000 0.7500000000
0.5000000000 0.0000000000 0.5000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
0.7500000000 0.2500000000 0.7500000000
0.7500000000 0.7500000000 0.2500000000
0.2500000000 0.2500000000 0.2500000000
!

for i in $(seq 400 50 700) #截断能从450-800,步数为50
do

#生成vasp静态计算的INCAR
cat > INCAR_static <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = .FALSE. (Projection operators: automatic)
ENCUT = $i (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)

Electronic Relaxation
ISMEAR = -5 (Gaussian smearing; metals:1)
#SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-06 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)

Ionic Relaxation
ISIF = 2 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V,4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
#KSPACING = 0.10
!

cp INCAR_static INCAR
echo "ENCUT = $i eV" ; time mpirun -np 16 vasp_std #vasp并行运行命令,根据系统自行修改

rm INCAR_static
rm INCAR
#提取计算得到的能量
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $E >>encut_energy.out
done

K点测试

是倒空间(动量空间)的基本构成点,只取在一个倒空间晶格向量的范围内来描述 k 。总能量计算是对布里渊区内的波函数进行积分,而积分是通过对部分特殊K点的求和完成的。因此,计算一个体系时,需要确定所用KPOINTS的大小。
k.sh如下:

{.line-numbers}
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#!/bin/bash
#生成vasp静态计算的INCAR
cat > INCAR_static <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = .FALSE. (Projection operators: automatic)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)

Electronic Relaxation
ISMEAR = -5 (Gaussian smearing; metals:1)
#SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-06 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)

Ionic Relaxation
ISIF = 2 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V,4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
#KSPACING = 0.10
!

cp INCAR_static INCAR

#产生计算所需POSCAR
cat > POSCAR<<!
Si8
1.0000000000
5.4687280000 0.0000000000 0.0000000000
0.0000000000 5.4687280000 0.0000000000
0.0000000000 0.0000000000 5.4687280000
Si
8
Direct
0.0000000000 0.0000000000 0.0000000000
0.2500000000 0.7500000000 0.7500000000
0.5000000000 0.0000000000 0.5000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
0.7500000000 0.2500000000 0.7500000000
0.7500000000 0.7500000000 0.2500000000
0.2500000000 0.2500000000 0.2500000000
!


for i in $(seq 6 3 15) #K点从6-24,步数为3
do

#生成KPOINTS文件,K点类型:Monkhorst-Pack (M) 或者 Gamma (G)
cat > KPOINTS <<!
Automatic mesh
0
Monkhorst-Pack
$i $i $i
0 0 0
!

echo "KPOINTS = $i" ; time mpirun -n 16 vasp_std #vasp并行运行命令,根据系统自行修改

rm KPOINTS
#提取计算得到的能量
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $E >>kpoints_energy.out
done

显示收敛情况的脚本

参考:https://blog.shishiruqi.com/2019/08/19/paili-energyforce/

  • 简单版:

grep F= OSZICAR |awk '{print $1,$5}'
grep F= OSZICAR |awk '{print $1,$5}' > conv.dat

  • 小脚本:
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    #!/usr/bin/env bash
    for i in *; do
    if [ -e $i/OUTCAR ]; then
    echo -e $i "\t" $(grep ' without' $i/OUTCAR |tail -n 1 | awk '{print $7}')
    fi
    done
  • 复杂:
    ediff.sh如下:(安装gnuplot)
    {.line-numbers}
    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
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    #!/bin/sh
    #lipai@mail.ustc.edu.cn

    begin=$1
    if [ $begin =='' ]; then begin=0; fi
    awk -v begin=$begin '/E0/{if ( i<begin ) i++;else print $0 }' OSZICAR >temp.e

    awk '/POSITION/,/drift/{
    if(NF==6) print $4,$5,$6;
    else if($1=="total") print $1 }' OUTCAR >temp.f

    awk '{if($4=="F"||$4=="T") print $4,$5,$6}' CONTCAR >temp.fix
    flag=`wc temp.fix|awk '{print $1}'`
    steps=`grep E0 OSZICAR |tail -1 |awk '{print $1}'`
    if [ flag != '0' ] ; then
    if [ -f temp.fixx ] ; then rm temp.fixx ; fi
    for i in `seq $steps`;do
    cat temp.fix >>temp.fixx
    echo >>temp.fixx
    done
    paste temp.f temp.fixx >temp.ff
    fi

    awk '{ if($1=="total") {print ++i,a;a=0}
    else {
    if($4=="F") x=0; else x=$1;
    if($5=="F") y=0; else y=$2;
    if($6=="F") z=0; else z=$3;
    force=sqrt(x^2+y^2+z^2);
    if(a<force) a=force} }' temp.ff >force.conv

    gnuplot <<EOF
    set term dumb
    set title 'Energy of each ionic step'
    set xlabel 'Ionic steps'
    set ylabel 'Energy(eV)'
    plot 'temp.e' u 1:5 w l t "Energy in eV"
    set title 'Max Force of each ionic step'
    set xlabel 'Ionic steps'
    set ylabel 'Force (eV/Angstrom)'
    plot 'force.conv' w l t "Force in eV/Angstrom"
    EOF

    if [ $steps -gt 8 ] ; then
    tail -5 force.conv >temp.fff
    tail -5 temp.e >temp.ee
    gnuplot <<EOF
    set term dumb
    set title 'Energy of each ionic step for the last few steps'
    set xlabel 'Ionic steps'
    set ylabel 'Energy(eV)'
    plot 'temp.ee' u 1:5 w l t "Energy in eV"
    set title 'Max Force of each ionic step for the last few steps'
    set xlabel 'Ionic steps'
    set ylabel 'Force (eV/Angstrom)'
    plot 'temp.fff' w l t "Force in eV/Angstrom"
    EOF
    rm temp.fff temp.ee
    fi

    rm temp.e temp.f temp.ff temp.fix temp.fixx

批量提交任务

适用于测试、有限位移法求力常数矩阵等情况。
简单: for i in *; do cd $i ; qsub sub4; cd $OLDPWD; done

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/sh
root_path='pwd'

#在每个job文件夹里进行静态计算。
for i in {001..1316} #按顺序从job-001116做静态vasp计算,根据生成的文件夹设定
do
cd job-$i
echo "job-$i", time
sbatch sub.sh
wait
cd ..
done

批量产生文件夹

i in {2..9}; do cp 0.01 0.0$i -r ; done ``` #批量复制0.01文件夹复制成0.02-0.09
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
```sed '3s/0.01/0.02/g' INCAR > INCAR_new```  #将INCAR中的第三行的0.01全部替换成0.02并输出位INCAR_new文件。
```sed -i '3s/0.01/0.02/g' INCAR``` # -i表示直接对源文件进行编辑,也就是说编辑之后源文件被新文件替换掉。
```for i in *; do sed -i "3s/0.05/$i/g" $i/INCAR ; done``` #将文件夹中所有的0.01替换成了文件夹序号。

我们主要对KPOINTS的文件的第四行进行批量操作,将1 1 1改成 2 2 2, 3 3 3 等。复制文件夹:
```for i in {1..6}; do cp ../ex03/0.01/ ${i}${i}${i} -r ; done```
显示行号:
```cat 333/KPOINTS -n ```
修改文件;
```for i in {1..6}; do sed -i '4s/1 1 1/$i $i $i/g' ${i}${i}${i}/KPOINTS ; done```

脚本:
```javascript {.line-numbers}
for i in 3RD.POSCAR.*
do
s=$(echo $i| cut -d "." -f 3) #代表.的3个字符
d=job-$s
mkdir $d
cp $i $d/POSCAR
cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d
pushd $d
popd
done

节约机时

grep User OUTCAR 后得出的结果被空格分成了4部分,时间信息在第4部分里面。User 或者 Elapsed都可以。

  • 一个合理的初始结构,可以避免很多意外的错误以及快速得到正确的结果。
  • 如果你设置的收敛标准比较严格,精度比较高,则每一个离子步需要更多的电子步数,需要的时间也会随着离子步的增加成线性关系增长。
  • ENCUT和KPOINTS文件需要测试,二者越大时间越久。
  • 这里我们就会需要一些其他的参数,NCORE和NPAR。
    NCORE:控制多少个核同时计算;
    NPAR:如何把计算任务分配到计算资源上面计算。
    它们之间的关系是:NCORE= 计算使用的核数 / NPAR。注意:这两个参数只能选取一个来使用:
    NCORE=单节点的核数,单节点的核数/2,单节点的核数/4…….个人使用经验是:NCORE = 单个节点核数 / 2 的时候,运行最省时间,设置也最方便。
    当然NCORE也需要测试。

消虚频

结构优化消虚频

X, Y, Z为体系中原子的坐标,也就是POSCAR中的内容。 dx, dy,dz为虚频对应的原子振动方向上的位移。我们将(dx, dy,dz)这个位移乘以一个介于0-1之间的校正因子,然后跟POSCAR中的坐标加起来即可。以z方向为例,校正因子设为0.1, 微调后的z方向坐标为: 10.709980 + (-0.676634)0.1。
脚本vib_correct.py如下:https://www.bigbrosci.com/2022/04/29/A29/

{.line-numbers}
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
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
from ase.io import read, write
import os

#获取虚频开始行
l_position = 0 #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
lines = f_in.readlines()
wave_num = 0.0
for num, line in enumerate(lines):
if 'f/i' in line:
wave_tem = float(line.rstrip().split()[6])
if wave_tem > wave_num: #获取最大的虚频
wave_num = wave_tem
l_position = num+2
# ASEPOSCAR
model = read('POSCAR')
model_positions = model.get_positions()
num_atoms = len(model)
#print(model_positions)

# 获取虚频对应的OUTCAR部分
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 72227248
#print(vib_lines)
vib_dis = []
for line in vib_lines:
#model_positions = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
vib_dis.append(vib_infor)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。
model.positions = new_positions

###保存结构
write('POSCAR_new', model, vasp5=True) # POSCAR_new是微调后的结构,用于下一步的计算(别忘了把POSCAR_new改成POSCAR)。

过渡态虚频

如果过渡态计算的时候出现好几个虚频,自己直接写个python脚本也能实现。不过使用ASE可以很方便地获取POSCAR中原子数目,保存原子固定的信息。参考如下。

{.line-numbers}
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
import os

#获取虚频开始行
l_position = 0 #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
lines = f_in.readlines()
wave_num = 0.0
for num, line in enumerate(lines):
if 'f/i' in line:
wave_tem = float(line.rstrip().split()[6])
if wave_tem > wave_num: #获取最大的虚频
wave_num = wave_tem
l_position = num+2
# 读POSCAR
with open('POSCAR') as f_pos:
lines_pos = f_pos.readlines()

# 获取虚频对应的OUTCAR部分
num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()])
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 72227248

model_positions = []
vib_dis = []
for line in vib_lines:
position = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
model_positions.append(position)
vib_dis.append(vib_infor)

# 微调结构
model_positions = np.array(model_positions)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。

###保存结构
f_out = open('POSCAR_new','w')
f_out.writelines(lines_pos[:8])
f_out.write('Cartesian\n')
for i in new_positions:
f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n')
f_out.close()

常用命令行

grep "33 81" -A3 FORCE_CONSTANTS | xargs echo -n| awk '{print $1,$2,$3+$7+$11}' >> trace.dat

计算原子间二阶力常数矩阵的迹

grep 'F=' OSZICAR|awk '{print $1, "\t"$3,"\t"$9}' >AIMD.dat

提取分子动力学温度能量数值

for i in 01 02 03 04 05;do cd $i && cp CONTCAR POSCAR && cd ..;done

CINEB计算重新计算

awk '{print $1,$2,$3+$4+$5,$6+$7+$8+$9+$10,$11}' PDOS > pdos.dat

pdos轨道加和

grep LOOP+ OUTCAR |awk '{sum+=$4} END {print "Average=",sum/NR}'

grep LOOP+ OUTCAR |awk '{sum+=$4} END {print "Sum=",sum}'

grep LOOP+ OUTCAR |wc -l

提取平均离子步时间,总时间,总步数

echo -e "4\n403\n49-96 121-144\n1\n"|vaspkit

vaspkit命令行

for i in *; do mkdir $i && cp INCAR KPOINTS POSCAR POTCAR sub.sh $i && sed -i "3s/0.05/$i/g" $i/INCAR ; done

批量修改INCAR

sed -i '5s/\<0\>/0.5/' filename
请将 filename 替换为您要修改的文件名。
5s/<0>/0.5/:
5:指定要修改的行号(第 5 行)。
s:表示替换操作。
<0>:匹配完整的数字 0,确保只替换单独的 0。
0.5:要替换成的新值。

  1. 插入的image只需要看一下有没有原子互相挨得很近,有的话调一下就行。
  2. 一开始计算先用IOPT=7,比IOPT=3快一点,然后LCLIMB=.FALSE.;算个一二十步之后再换IOPT=1或2,LCLIMB=.TRUE.
  3. NEB计算不需要收敛,从经验来看也极少能碰到收敛的情况。主要是找到切线力尽可能等于0的点,然后用dimer算,基本都能找到过渡态。
  4. nebresult.pl没用过,看结果“forces and energy:”下面的数据的第四列应该就是切线力,正常的计算切线力是由正值变到负值,或者由负值变到正值,然后找接近于0的点。你图里的数据只有正值,因此要么初末态构型有问题,要么你要找的点在0和1两个image之间。
    不过1号image的切线力已经很接近0了,可以尝试直接做过dimer
  5. 跑个一二十步之后停止计算,‘mv CONTCAR POSCAR’之后修改INCAR再继续算。脚本都是线性插值,因此初始构型一般都很不合理,上来就用IOPT=2的话收敛慢,甚至可能出错。跑个几十步等构型相对合理一些之后,再换IOPT=2。
  6. dimer是vasp里面自带的算法,详细用法手册里有。

计算步骤

  1. 初末态结构优化
    建立两个文件夹ini(初态的意思),fin(末态),每个文件夹放入vasp计算必备的四个文件(INCAR,POSCAR,KPOINTS,POTCAR),其中的两个POSCAR对应未优化的初末态。确保两个文件夹里面除POSCAR外,其他文件完全一样。

对于两个POSCAR,注意每一行的原子一一对应。若有固定位置的原子(比方做表面计算,要固定底部1-2层原子),最好让两个POSCAR里面的固定原子位置完全相同。
用dist.pl检查两个优化后结构(两个CONTCAR)的相似程度(每个对应原子的初末态距离的平方和,再开根号)。

具体命令为: dist.pl ini/CONTCAR fin/CONTCAR
若返回值小于5A,一般可以进行下一步。
2. 插点
插点的数目取决于前面dist.pl的返回值,一般插点数目可取(dist.pl返回值/0.8)。

具体插点的命令: nebmake.pl ini/CONTCAR fin/CONTCAR N

最后的N表示插点数目。插入点的算法为线性插值,详情请进前面给的vtst脚本链接。我这里设置N为3,执行完命令后出现文件夹00,01,02,03,04.五个文件夹里面都有nebmake.pl产生的POSCAR。

00表示初态,里面放的是ini/CONTCAR,04表示末态,放的是fin/CONTCAR, 01、02、03是插入的点。

输出中提示讲初末态对应的OUTCAR复制到对应的文件夹中,以便后续数据分析。
检查插入点的合理性:

输入命令 nebmovie.pl 0

参数0表示用POSCAR生成xyz文件;还可取1,为用CONTCAR生成。

可以看到每个文件夹内多了一个后缀名为xyz的文件。将这几个xyz文件合成一个movie.xyz文件即可在VMD中动画演示。我更倾向于使用jmol(或者ovito)来查看xyz文件,不过VMD可能大家更熟悉。

用nebmovie.pl,没有直接生成movie.xyz是因为从官方主页下载的脚本默认当使用nebmovie.pl后面不带参数或者参数为0的时候(即使用POSCAR产生xyz),不输出movie.xyz。 你想让脚本自动生成movie.xyz而不是自己去合并各个文件夹里面的xyz文件,需要自行修改nebmovie.pl文件。很简单,把倒数第二个if语句整个用#注释掉或者直接删掉即可。

也可以直接用VESTA查看01的POSCAR。

进一步可用命令nebavoid.pl 1,确保中间插入的点每一个原子间距都大于1A。该命令的参数1表示最小允许间距,可取小数。
原子间距太小说明结构有问题。要么是因为初末态POSCAR里的原子没有对应好,要么是因为nebmake.pl线性插值的方法不适用于估计反应路径路径。如果是前者,调整好对应原子,也就是CONTCAR里面相应行进行调换;如果是后者,用Materials studio之类的软件调整好结构,再放入对应文件夹。
3. 其它文件的准备
在当前目录下面放入KPOINTS,POTCAR及INCAR文件。要求KPOINTS和POTCAR与ini、fin文件夹中的一样;INCAR中的基本参数也与初末态计算的INCAR一样,另外再加入NEB计算的相关参数。
INCAR再多说两句。优化算法若想使用vtst的,需要如上设置IBRION=3,POTIM=0,IOPT设置见本文第四个图。若使用vasp自带的优化器,可以使用IBRION=1或者3,不要取2;POTIM取合适的值(0.01-0.5范围内去尝试)。

这样就可以提交NEB任务了。
4. 后处理
计算过程中,可以用nebef.pl命令(不带参数)查看计算收敛情况。
输出中,第二列即为最大受力(force of images in the neb),第三列为相应结构的能量。

前面INCAR中EDIFFG设置为-0.02(一般结构优化可取-0.02,NEB计算取-0.03),表示当所有结构最大力小于0.02eV/A时,结束计算。

另一个可以用来观察收敛情况的命令是nebbarrier.pl(同样不带参数)。该命令没有输出到屏幕的内容,而是生成neb.dat文件。
neb.dat文件第二列表示距离(即临近两结构的dist.pl的计算结果),第三列表示能量(以初态能量为参考值),第四列为力(forces along the neb)。

EDIFFG参数对应的力是nebef.pl输出中的force of images in the neb.

计算完成后,使用命令nebresult.pl.

nebresult.pl做的事情如其所输出表明的,完成了nebbarrier.pl, nebspline.pl, nebef.pl, nebmovie.pl, nebjmovie.pl, nebconverge.pl还有对各文件夹中的OUTCAR打包压缩。生成了很多文件。其中mep.eps是以dist.pl距离为横坐标,能量为纵坐标画出的能势垒图,如下:
该文件一般查看图片的软件打不开,推荐使用EPS/PS viewer. https://epsviewer.org/download.aspx

生成的spline.dat文件是对上面几个点的拟合曲线数据,已经有了上图,这些数据就意义不大了。

生成的vaspgr文件夹内是各个插点结构的收敛图(同样可用EPS/PS viewer打开),如下vaspout1.eps是01结构收敛信息:

加压

可通过施加静水压强、调整晶格常数的方法实现结构压缩

vi INCAR 增加标签 PSTRESS=50
(单位为 0.1 GPa)

其他不变

加电场

在vasp中根据自己所需要添加电场的场强更改INCAR中的如下参数:

1
2
3
LDIPOL=.TRUE.
IDIPOL=1/2/3/4
EFIELD= ?

对于IDIPOL参数的设置,需注意:对于IDIPOL=1-3,偶极矩将分别平行于第一、第二或第三晶格矢量的方向计算。总能量的修正计算为当前超单元中单极子/偶极子和四极子与放置在超单元中的同一偶极子之间的能量差,对应的晶格矢量接近无穷大。此标志应用于板计算。对于IDIPOL=4,将计算所有方向上的完整偶极矩,并将总能量的修正计算为当前超级单体中的单极/偶极/四极和放置在真空中的同一单极/偶极/四极之间的能量差,使用该标志计算孤立分子。

此外,EFIELD 的设置需要注意单位的换算。eV/Å(电子伏特/埃)计算完毕之后通过后处理软件分析所关注的性质即可。

状态方程拟合

我们需要对结构进行优化才可以获取稳定的晶格参数信息。有两个方法可以实现:

1 Birch-Murnaghan状态方程拟合,
2 VASP计算中通过调节ISIF参数直接优化Bulk。
直接优化前面已经介绍,下面我们先讨论一下第一个方法:BM方程拟合。

参考:大师兄科研网:https://www.bigbrosci.com/2018/02/02/ex33/

科学网耿华运的博文:https://blog.sciencenet.cn/blog-3469498-1280203.html

科学网赵志强的博文,关于拟合:https://wap.sciencenet.cn/blog-3388193-1152915.html

vaspkit:http://vaspkit.cn/index.php/48.html

大致流程(以大师兄提供脚本为例)

  1. 获得一定缩放系数的结构:for i in 0.95 0.96 0.97 0.98 0.99 1.01 1.02 1.03 1.04 1.05; do cp 1.00 $i ; sed -i "2s/1.0/$i/g" $i/POSCAR ; done

  2. 分别进行结构优化:for i in *; do cd $i; 提交任务命令; cd $OLDPWD; done

  3. 提取能量-体积以及能量-压力、压力-体积的数据:for i in *; do echo -e $i "\t" $(grep ' without' $i/OUTCAR | tail -n 1| awk '{print $7}'); done > data

  4. 拟合能量体积曲线,获得平衡晶胞参数

  • 这个流程可通过批量计算脚本实现,参考上面链接。另外还可以获得能量随压力的变化来模拟相变。

分数占据

生成disorder或者掺杂构型的软件汇总:

https://mp.weixin.qq.com/s/vYar9UcYWAZgazBTgqUGEA
https://github.com/jichunlian/disorder
supercell: https://orex.github.io/supercell/
ATAT: https://www.brown.edu/Departments/Engineering/Labs/avdw/atat/
vaspkit:https://vaspkit.com/

disorder

安装

1
2
3
4
5
6
7
8
$ unzip disorder-main
$ cd disorder-main
$ make disorder (Only disorder is compilied)
$ make supercell (Only supercell is compilied)
$ make or make all (Both disorder and supercell are compilied)
$ make clean (rm *.o *.mod disorder supercell)
export PATH=$PATH:disorder_installation_path/disorder-main/bin
(Add this to the ~/.bashrc file as you wish)

Note: The default compiler is ifort, and you can modify the Makefile
in the src directory to use another compiler (e.g. gfortran).

需要准备两个文件SPOSCAR和INDSOD,以examples的1和2为例:

例子1:SnxPb1-xTe

  1. 扩建超胞,保证整除
    1
    2
    3
    4
    5
    6
    7
    8
    9
    Te Pb
    1.0
    13.1323995589999996 0.0000000000000000 0.0000000000000000
    0.0000000000000000 13.1323995589999996 0.0000000000000000
    0.0000000000000000 0.0000000000000000 13.1323995589999996
    Te Pb
    32 32
    Direct
    * * *
  2. 准备INDSOD文件
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    &input
    nsub = 2
    subs = 16,16
    symb = 'Sn','Pb' ! The quotes is unnecessary for the ifort compiler
    site = 2
    prec = 1D-5
    ! lpro = .true.
    ! lpos = .true.
    ! leqa = .true.
    ! lspg = .true.
    lcfg = .false.
    /
    官方解释:
    1
    2
    3
    4
    5
    6
    7
    nsub-tag:nsub = 2 ~ 5; Default: nsub = 2.该位置取代原子种类。本例中是2种。
    subs-tag:subs = integer, integer, …; Default: None。对 应的原子数目。总和加起来应该等于该位点在SPOSCAR中的数目。
    symb-tag:symb = character, character, …; Default: None。元素化学符号。空位用Kw,即汉语kong wei。本例是Sn Pb单引号。
    site-tag:site = integer; Default: 1。所取代的位置在SPOSCAR中的位置。本例是第2.
    prec-tag:prec = real; Default: 1D-5。当软件产生空间群与其他软件识别不一致时调整该参数,一般不用修改。
    lpro-tag:lpro = logical; Default: .false.是否消除重复参数。对于大体系可以打开。
    lpos、leqa、lspg、lcfg:分别控制输出文件POSCAR-X、EQAMAT、SPGMAT、CFGMAT。
  3. 运行disorder。输入disorder即可。
  4. 输出文件:

ATAT软件

简介

  • 什么是ATAT?

其全称为 Alloy Theoretic Automated Toolkit (ATAT)
引用官网的一句话,ATAT is a generic name that refers to a collection of alloy theory tools developped by Axel van de Walle, in collaboration with various research groups and with various sources of financial support.

  • 如何安装ATAT:

http://brown.edu/Departments/Engineering/Labs/avdw//atat/atat3_36.tar.gz下载ATAT源码,解压后,make all开始编译;之后make install进行程序和脚本定向安装;之后在home中打开./bashrc添加ATAT安装路径于环境变量当中(添加 export PATH=$PATH:/xxxx),并source ./bashrc即可。

进入安装包的glue目录,里面有与各种dft软件的算例的接口设置。例如进入vasp可以看到它的接口设置脚本ezvasp,由于make install已经将它拷贝到了安装目录中,所以在bin中打开ezvasp,并设置相关可执行程序,包括赝势文件的名称和路径即可。

配置~/.ezvasp.rc

1
2
3
#!/bin/csh
set VASPCMD="/path/vasp_std# 设置vasp可执行程序
set POTLDA/PBE="/path/to/pot_lad" # 设置为赝势实际的路径
  • ATAT能做啥?

ATAT 是Walle及其合作者开发的一个利用第一性原理结合MonteCarlo模拟计算相图及热力学性能的软件包。该软件包首先采用团簇展开(cluster expansion),用第一原理方法计算一系列构型的形成能,然后用这些能量作为Monte Carlo模拟的哈密顿量,计算材料的热力学性能和高熵合金的mcsqs建模。

MAPs

软件运行

一)运行MAPS:

  1. 准备lat.in文件

2)maps -d & #maps命令运行,等待后面的指令

3)touch ready

#等待十几秒钟,出现“0”文件夹,是合金中比例为0的纯相结构。

4)此时需要利用第一性原理软件计算该结构的能量

#(参见下部分和VASP软件的交互,可无需手动输入),现在手动输入能量,理解一下计算过程,假设该结构的能量已经知道,在0这个文件夹下,echo 1.1 > energy, 并删除wait。继续touch ready,就会生成1这个文件夹,里面是比例为1的另一个纯相结构,echo 2.2 > energy, 删除wait。touch ready,生成2这个文件夹,一直这样循环下去,直至达到相应的精度,最终停止输入touch stop。

二)运行MAPS与VASP交互:

1)需要lat.in,vasp.wrap文件 ( #vasp的INCAR参数)

2)其余与上面运行MAPS步骤相同,直到第4步骤计算能量,进入该文件夹,运行runstruct_vasp。

大概率会报错,原因有几个:

a. 常见为~/.ezvasp.rc文件配置问题,赝势位置,vasp运行路径出错。

可以看运行命令后,是否有POTCAR文件产生,如果没有,基本就是赝势位置问题,如果POTCAR文件正常,就是vasp没有运行。

b. vasp不运行,在集群上面运行的软件,ATAT软件给了集群上面的运行办法, pollmach命令,将runstruct_vasp命令写在脚本里是可以的。

最终正确版本运行步骤:

1)准备lat.in,vasp.wrap文件,脚本文件。

lat.in文件

vasp.wrap文件

sub.sh文件

2)提交脚本文件即可

说明:

  • 缺少.machines.rc文件的提示可忽略,该文件用于跨节点的并行
  • maps不断产生结构和构建团簇展开模型,-d表示使用默认参数,直接输入maps可查看可用参数
  • runstruct_vasp调用vasp计算所产生的结构的能量,每个目录包含一个结构及其计算结果,每个目录下的energy文件即从vasp的计算中读取的能量
  • maps.log中记载了现有的团簇展开模型的状态,包括对基态构型的预测和LOOCV等

什么时候停止计算:

  • 团簇展开模型能够成功预测基态 (查看maps.log)
  • LOOCV足够小 (<0.025 eV suggested by ATAT manual)
  • 有效团簇相互作用随着团簇距离的增大而衰减

以上的三个判据并不是绝对的,判断一个团簇展开模型是合理的还可以有其它方式。

停止计算:

1
2
$ touch stop
$ touch stoppoll

重启计算: 与正常开始一个计算相同,如果之前的计算不是正常结束则需要删除一些文件

1
2
3
rm maps_is_running pollmach_is_running ready
maps -d &
pollmach runstruct_vasp

常用命令

  • corrdump -2= -3= -4= 产生团簇

  • getclus 查看团簇

    • 第一列为团簇的阶数,第二列为团簇直径,第三列为团簇的多重度
    • clusters.out中有团簇格点位置的详细信息
  • clusterexpand -e -cv energy 得到团簇展开系数

    • 此时产生的ECI文件为energy.eci,而maps产生的ECI为eci.out,后者是形成能的团簇展开结果,前者是直接对vasp能量(即energy文件)的团簇展开
    • 如果需要对其它性质做团簇展开,例如带隙,则可以在每个结构的目录下创建一个bandgap文件,其内容即为带隙的大小,然后运行clusterexpand -e -cv -pa bandgap以得到带隙团簇展开的系数bandgap.eci,其中-pa表示所展开的量已经是一个平均的量

检查计算结果

  • VASP是否正常结束
    • ATAT会自动检查一些错误(checkerr_vasp)
    • 但如果结构优化的离子步数不够大,不会有提示,因此建议在vasp.wrap中将NSW设置大一些
  • 晶胞的应变是否很大
    • checkrelax (<0.1 suggested by ATAT)

结果可视化

  • mapsrep 会调用gnuplot作图并产生相应的脚本mapsrep.gnu

输出文件

  1. maps.log:开始可能的warnings:
    ’Not enough known energies to fit CE’
    ’True ground states not = fitted ground states’
    ’New ground states predicted, see predstr.out’
    随着结构增多:
    ’Among structures of known energy, true and predicted ground states agree.’
    ’No other ground states of xx atoms/unit cell or less exist.’
    最后一行有crossvalidation score (in eV/atom),该值要小于0.025。精确小于0.001.

  2. fit.out:拟合的结构,一个结构一行,共6列,分别代表concentration、 energy 、fitted_energy (energy-fitted_energy)、 weight 、index

  3. predstr.out:预测的能量,一个结构一行,共5列,分别代表concentration、 energy 、predicted_energy、index、status.其中status包括status is either
    b for busy (being calculated),
    e for error,
    u for unknown (not yet calculated) or
    g if that structure is predicted to be a ground state (g can be combined with the above).

  4. gs.out:基态能量,一个结构一行,四列,分别是concentration、 energy、 fitted_energy、 index

  5. gs_str.out:基态的结构,格式和str.out一样。每一个结构后面会有end和一个空行隔开。

  6. eci.out:eci值对应的cluster在clusters.out文件。

  7. clusters.out:对于每一个团簇,第一行是多重度 multiplicity,第二行是团簇维度cluster diameter,第三行是团簇的格点数,下一行是格点对应的坐标。最后空行隔开每一个团簇。

  8. str.out:在每个数字文件夹下面,和晶格输入文件相似,不同是坐标写成3*3矩阵,一个位置只写一个原子。

  9. ref_energy.out:用于计算形成能的参考能量。

后处理

一系列用于分析ATAT计算结果的python脚本,GitHub地址为:https://github.com/xu-xi/clusterexpansion可以通过 git clone 下载到本地服务器。用“-h”命令可以查看每个脚本的帮助文档。最常用的几个脚本是:

1)cv_scan.py,测试CE模型对团簇截断半径的收敛性。

2)lc.py,测试CE模型对训练集大小的收敛性。通常~80个构型就可以得到比较精确的团簇展开模型。

3)ce_plot.py,这个脚本绘制的是拟合能量-计算能量分布图,默认-p=energy(调用的是每个数字目录下的energy文件)。energy文件是ATAT完成一个构型的DFT (relax+static) 计算后自动生成的能量文件。在实际应用中,拟合能量-计算能量分布图需要的往往不是energy而是formation energy,所以就需要编写脚本手动在每一个数字目录下生成一个形成能文件formation_energy,该值需要将目标目录下energy的值与0文件夹的energy和1文件夹的energy作线性变换得到,再写入到每一个数字目录中,读者可自行尝试。下面的脚本仅供参考,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#!/usr/bin/env python
import os,subprocess,numpy,argparse,sys,subprocess,math
energy_0 = float(subprocess.check_output('cat ./0/energy', shell = True))
energy_1 = float(subprocess.check_output('cat ./1/energy', shell = True))
for item in os.listdir(os.getcwd()):
fullpath = os.path.join(os.getcwd(),item)
if os.path.isdir(fullpath):
os.chdir(fullpath)
if not os.path.isfile('error'):
sc_atom_number = len(file('str.out').readlines())-6
# print sc_atom_number
energy = float(subprocess.check_output('cat energy', shell = True))
energy /= sc_atom_number
#energy *= 5 normalize to unit cell,5 atoms in CsPbI3
sc_x = float(subprocess.check_output('corrdump -pc -l=../lat.in',shell=True).strip().split()[0])
formation_energy = energy - (1 - sc_x) * energy_0 - sc_x * energy_1
# print formation_energy
with open('formation energy','w') as f:
f.write('%.7f\n' % formation_energy)
os.chdir('../')

注:这里如果报错:NameError: name ‘file’ is not defined. Did you mean: ‘filter’?将file改成open

得到formation_energy文件后,执行命令ce_plot.py -p=formation_energy -s=”formation energy”,还可以用Origin手动绘制出形成能-掺杂浓度分布图;在得到了比较精确的CE模型后,可以进行蒙特卡罗模拟,得到团簇关联函数、系综平均能量和方差随温度变化示意图。

热力学计算

必要的输入文件:

  • lat.in
  • clusters.out
  • eci.out
  • gs_str.out
emc2使用示例

蒙特卡洛模拟

1
emc2 -T0=300 -T1=3000 -dT=100 -cm -x=0 -keV -gs=-1 -er=30 -aq=0 -dx=1e-5 -sigdig=12
  • -cm表示正则系综
  • -x其实是点团簇函数,若c表示格点的占据浓度,则x=2c-1
  • -sigdig=12是为了表示有足够的位数输出能量的方差以判断相变
  • -aq-dx是自动确定MC平衡步数和取样步数的算法,其原理可参考文献Ref.1,也可以用-eq-n来直接指定步数,需要做收敛性测试
  • 主要输出文件为mc.out,输入emc2 -h可查看每一列所表示的内容
  • 输出文件mcsnapshot.out为最后一个模拟温度下的MC快照。如果想得到某个温度下的snapshot,可以只设置T0而不设置T1和dT。
    1
    2
    3
    $ cp mcsnapshot.out str.out
    $ runstruct_vasp -nr
    $ vesta POSCAR

相图计算

基本原理可参考文献Ref.1,Martin Baker的笔记提供了更多的例子和细节

例子:

1
phb -dT=10 -ltep=1e-3 -er=30 -gs1=0 -gs2=1 -o=phb.out -keV -dx=1e-5

构建SQS/SQoS

SQS/SQoS是对无序体系的代表性结构,其中SQS假定原子完全随机的占据,而SQoS是SQS的扩展,可以考虑一定的短程有序性,其基本原理可参考相关文献Ref.2

必要文件:

  • rndstr.in
  • clusters.out

rndstr.in的一个例子:

1
2
3
4
5
6
7
8
9
10
11
4.1000 0.000000 0.000000
0.000000 4.1000 0.000000
0.000000 0.000000 4.1000
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
0.000000 0.500000 0.500000 O=0.6667,N=0.3333
0.500000 0.000000 0.500000 O=0.6667,N=0.3333
0.500000 0.500000 0.000000 O=0.6667,N=0.3333
0.000000 0.000000 0.000000 Ba
0.500000 0.500000 0.500000 Ta

产生团簇文件clusters.out

1
corrdump -l=rndstr.in -clus -noe -nop -2= -3= -4=

说明:

  • -noe-nop表示忽略空团簇和点团簇
  • getclus查看所产生的团簇。如果团簇数目过少,则后续产生SQS有可能很容易完全匹配团簇函数,导致SQS实际上不够随机;而如果团簇过多,则有可能导致近邻的团簇匹配得不够好,而通常近邻的团簇是更为重要的,可以通过-wd加上权重以避免这种情况

运行mcsqs -n= &
说明:

  • -n设定SQS晶胞的原子数
  • 可通过-tcf设定目标团簇函数(可从mc.out中提取),即SQoS方法
  • mcsqs不会自动停下来,可用touch stopsqs停止
  • 输出文件mcsqs.log中有对团簇函数匹配的优化记录
  • 输出文件bestsqs.out为目前找到的最佳的SQS,可将其复制为str.out,然后运行runstruct_vasp -nr生成vasp的输入文件用于后续计算

MCSQS

一、原理

无序态的关联函数已知,因此只需要备选态的关联函数和无序态关联函数差别最小就可以,即让下方公式中的目标函数最小。

二、运行步骤

step1. 准备rndstr.in文件

1 1 1 90 90 90

1 0 0

0 1 0

0 0 1

.0 .0 .0 Cr

.0 .5 .5 Ni=0.333333,Fe=.333333, Co=.333334

.5 .0 .5 Ni=0.333333,Fe=.333333, Co=.333334

.5 .5 .0 Ni=0.333333,Fe=.333333, Co=.333334

#第一行也可以改成3*3矩阵,类似于POSCAR中晶格常数格式。这里一般不用实际的晶格常数,用1替代。

用1对后面的参数设置更容易。

step2. 运行命令
corrdump -l=rndstr.in -ro -noe -nop -clus -2=1.1
该命令通过读取rndstr.in文件中的结构,进而计算结构的对称性和团簇等信息,输出clusters.out文件。

#需要注意的是‘-2’这个参数,是用于计算关联函数的截断距离。如对于FCC构型而言,晶格常数为1,第一最近距离为2^0.5/2=0.71,第二距离为1,第三距离为1.5^0.5=1.2。所以这里设置阶段距离在1和1.2之间。一般考虑‘-2’参数即可。(注:更高的-2,-3参数会产生更无序的结构,但是需要更长的时间,一般晶格常数为1,这里选取1.1就可以)

step3. 运行命令 mcsqs -n XX
该命令主要输出sqscell.out文件。XX表示体系的原子个数,该数值一定是上述结构原子的倍数,否则会出错。

#比如上述构型,自身CrNi有4个原子,其中Ni元素与Fe,Co构成合金,因此XX至少为12个原子。该步骤运行一下,就可以进行停止。输出文件需要关注sqscell.out文件,该文件给出了XX个原子扩胞后形成的不同晶格构型。可通过挑选,删除一些不太可能的结构,比如x,y,z三个方向长度差别过大,形成的角度较小等对VASP或后续计算不友好的结构。之后用新的sqscell.out文件运行第四步骤。

step4. 运行命令mcsqs -rc
输出bestsqs.out文件和bestcorr.out文件。其中bestsqs.out就是无序结构。

#通常而言没有最完美的无序结构,这个命令会一直运行下去,可以肯定的是,再一次生成的结构比上一次生成的结构要优异,所以什么时候停止取决于自身研究的精度。

层间结合能、结合力

贴一篇文章:https://news.sciencenet.cn/htmlpaper/2022/4/202242914585844872440.shtm

Graphite interlayer distance

vasp上的教程:

准备graphite结构,

INCAR:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 20           
LVDW_EWALD =.TRUE.
NSW = 1
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

提交脚本:
To run this example, execute the run.sh bash-script:

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

rm results.dat

for d in 6.5 6.6 6.65 6.7 6.75 6.8 6.9 7.0
do
cat>POSCAR<<!
graphite
1.0
1.22800000 -2.12695839 0.00000000
1.22800000 2.12695839 0.00000000
0.00000000 0.00000000 $d
4
direct
0.00000000 0.00000000 0.25000000
0.00000000 0.00000000 0.75000000
0.33333333 0.66666667 0.25000000
0.66666667 0.33333333 0.75000000
!
$vasp_std
cp OUTCAR OUTCAR.$d
energy=$(grep "free ene" OUTCAR.$d|awk '{print $5}')
echo $d $energy >> results.dat
done

The optimal length of the lattice vector c normal to the stacking direction is determined in a series of single point calculations with varied value of c (all other degrees of freedom are fixed at their experimental values).

The computed c vs. energy dependence is written in the file results.dat and can be visualized e.g. using xmgrace. The optimal value can be obtained using the attached utility (python with numpy or Numeric is needed):

./utilities/fit.py results.dat
This should yield:

1
2
3
4
5
6
7
200 iterations performed
Ch-square: 4.30305519481e-09
---------

E0(eV): -37.433456779
d0(A): 6.65603352689
The computed value of 6.66 Å agrees well with experiment (6.71 Å).

Graphite MBD binding energy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 202           
LVDWEXPANSION =.TRUE.
NSW = 1
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

KPOINTS

Graphite:

1
2
3
4
5
Monkhorst Pack
0
gamma
16 16 8
0 0 0

Graphene:

1
2
3
4
5
Monkhorst Pack
0
gamma
16 16 1
0 0 0

Running this example
To run this example, execute the run.sh bash-script:

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

# Here the work starts
rm results.dat

drct=$(pwd)

for i in graphene graphite
do
cd $drct/$i
$vasp_std
done

cd $drct

# obtain total energy for graphite
en2=$(grep "free ene" graphite/OUTCAR |tail -1|awk '{print $5}')

# obtain total energy for graphene
en1=$(grep "free ene" graphene/OUTCAR |tail -1|awk '{print $5}')

# compute interlayer binding energy (eV/atom)
deltaE=$(echo print $en2/4 - $en1/2 |python)

echo "Binding energy (eV/atom): " $deltaE >results.dat

Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)

The computed value of 0.050 eV/A is now fairly close to the RPA reference of 0.048 eV/atom [1].

Graphite TS binding energy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 20            
LVDW_EWALD =.TRUE.
NSW = 1
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

To run this example, execute the run.sh bash-script:

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

rm results.dat

drct=$(pwd)

for i in graphene graphite
do
cd $drct/$i
$vasp_std
done

cd $drct

# obtain total energy for graphite
en2=$(grep "free ene" graphite/OUTCAR |tail -1|awk '{print $5}')

# obtain total energy for graphene
en1=$(grep "free ene" graphene/OUTCAR |tail -1|awk '{print $5}')

# compute interlayer binding energy (eV/atom)
deltaE=$(echo print $en2/4 - $en1/2 |python)

echo "Binding energy (eV/atom): " $deltaE >results.dat

Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)

Even though the TS method predicts a reasonable geometry (see the Graphite interlayer distance example) it overestimates the energetics strongly: the computed binding energy of -0.083 eV/atom is too large compared to the RPA reference of 0.048 eV/atom. This overestimation is - at least in part - due to neglecting the many-body interactions (see example Graphite MBD binding energy).

结构优化完成后,我们已经可以从输出文件得到很多性质。这里我们先关注能量。

VASP输出的能量是什么

在优化完的文件夹下输入

{.line-numbers}
1
2
[scf1355@ln151%nc-l 2h]$ tail -1 OSZICAR 
3 F= -.44829871E+02 E0= -.44829871E+02 d E =-.169526E-03 mag= -0.0000

这里的F和下面free energy含义相同,E0代表体系内能。
另外在OUTCAR中还可以找到这些
{.line-numbers}
1
2
3
free  energy   TOTEN  =       -44.82987077 eV
energy without entropy= -44.82987077 energy(sigma->0) = -44.82987077
atomic energy EATOM = 1542.31592482

free energy指的是亥姆霍兹自由能(Helmholtz free energy),在VASP手册中被称为自由电子能(free electronic energy),它可以用F表示: F = E-TS+PV .而F和E对应上面的F和E0。
TOTEN和without entropy分别代表考虑和未考虑熵的贡献时自由能的值。而energy(sigma->0)为展宽趋近于零时的自由能,可以外推得到: energy(sigma->0) = E0 =1/2(F+E)

其中上面的展宽对应于INCAR中设置的ISMEAR参数,详情可以查看手册或者参考 https://zhuanlan.zhihu.com/p/384626631。
https://www.bilibili.com/video/BV19L4y1i7mC/?vd_source=17084cc867ff3bb86398ff1a79b443f2
简单理解费米面内电子全占据,费米面外没有电子。这种阶梯函数就使得,随着K点数目变化,K点处积分值变化很大,不易收敛。于是我们需要一个连续函数,代替该阶梯函数进行积分,这种替代方法就叫smear。在VASP计算中,只有在ISMEAR=-5时,才能严格满足费米能级的定义;而在其他情况下,允许电子在费米能级附近的宽度内“波动”,即将阶梯函数替代为一个平滑的函数,从而在不破坏和的精度的情况下获得更快的收敛速度。

明白非常重要的一点,从OUTCAR中直接提取的能量在不同赝势下不具备重复性,是不能够直接写到论文中去的。一般计算相对能量的时候只要保持取一样的值即可 比较绝对能量没啥意义。

形成能

形成能很简单,参考高中化学形成焓的计算可以知道对于反应A+B=C+D,那么形成能=E(C)+E(D)-E(A)-E(B)
同样对于MoS2的生成焓我们可以知道: $\Delta$ H=E(MoS2)-E(Mo)-2E(S)
因此还需要计算Mo单质和S单质的单原子能量。计算过程和1.1完全一致,只需要获得结构,利用vaspkit生成KPOINTS和POTCAR文件,INCAR可以不用改。能量用F或者E0或者上面说的哪一项都可以,但要保持一致。代入计算即可。

除此之外还可以计算产生一个Mo的空位的形成能: $\Delta$ H=E(空位,MoS2)+E(Mo)-E(MoS2)
但是此时没有考虑缺陷的电荷。
可以计算分解能,剥离能,滑移能,表面能,偏聚能,吸附能等等。但是具体结构的优化会有一些处理上的不同。

单点能SP

参考: https://www.bilibili.com/read/cv10064345

单点计算是指对结果不做优化,直接输出能量。因此优化只会进行一步。

设置NSW = 0 也可以不设置,因为默认值就是0。 其他文件和参数不需要改变。提交即可。

对于一个结构在结构优化(NSW不等于0)时,第一步输出的1F和单点能的结果是一样的。

自洽计算

自洽计算其实也是单点计算,静态自洽计算前需要提供 较稳定体系 的晶格结构信息(即结构优化完的CONTCAR),从而通过电子自洽计算,通过自洽迭代求解薛定谔方程(微观中描述电子的状态,相当于宏观的运动方程)) 完整地计算出体系基态下费米能级(准确值,The fermi level location is accurate only in self-consistent calculation.)、 电子的波函数、电荷密度等信息,可以直接分析原子间的键合作用,也可以在非自洽之后进一步分析晶体的电子结构和材料的相关性质。

因此对于电子结构计算之前必须进行一步自洽计算优化电子分布,以便加速后续计算。

{.line-numbers}
1
2
3
NSW=0  不再进行原子迟豫
LWAVE=.TRUE (需要在WAVECAR中输出波函数)
LCHARGE=.TRUE (需要在CHGCHGCAR中输出电荷密度)

气体优化

在涉及单质结构优化时,可能会涉及到气体分子的优化。

  • 一般情况下选择建一个相对大的格子,里面放入一个气体分子,格子的大小保证分子间不会影响。其他的和前面结构优化保持一致即可。
  • 另一个办法就是找数据库, ChemScpider。网址:http://www.chemspider.com

IVDW参数

有两大类的范德华斯作用修正:

  1. 基于半经验的,包括了D2, D3, D3-BJ, TS, TS+SCS等等,这些都是在常用的交换关联泛函比如PBE的基础上,考虑色散力的作用,在PBE计算出来的总能基础上增加了额外半经验项,这一项需要设置一些参数,至于是哪种半经验公式,由IVDW的设置来决定,具体的参考VASP的手册。
    范德华力计算方法,在DFT能量计算基础上增加范德华力修正。IVDW=10:DFT-D2方法;IVDW=11:DFT-D3方法。推荐首选更新的DFT-D3方法。
  2. 另一类就是vdW-DF,也就是修改交换关联项来考虑范德华修正。包括你所提到的optPBE-vdW, optB88-vdW, and optB86b-vdW等。它们的设置由修改GGA赋值,以及额外相关的一些参数来设置。 这些在vasp手册上有明确的说明。
    乱谈DFT-D:http://sobereva.com/83
    DFT-D色散校正的使用:http://sobereva.com/210
    大体系弱相互作用计算的解决之道:http://sobereva.com/214
    LUSE_VDW:https://blog.csdn.net/Yaris_liu/article/details/124683272 https://zhuanlan.zhihu.com/p/47619708
    包含vdW-DF,optPBE-vdW, optB88-vdW, optB88b-vdW, vdw-DF2。需要在所在的文件夹下放置一个名叫vdw_kernel.bindat的文件。
  • vdW-DF:采用该范德瓦尔斯修正需在INCAR中加上以下开关:
    {.line-numbers}
    1
    2
    3
    GGA=RE
    LUSE_VDW=.TRUE.
    AGGAC=0.0000
  • optB88-vdW:采用该范德瓦尔斯修正需在INCAR中加上以下开关:
{.line-numbers} GGA
1
2
3
4
PARAM1=0.1833333333
PARAM2=0.2200000000
LUSE_VDW=.TRUE.
AGGAC=0.0000
  • optB86b-vdW:采用该范德瓦尔斯修正需在INCAR中加上以下开关:
{.line-numbers} GGA
1
2
3
4
PARAM1=0.1234
PARAM2=1.0000
LUSE_VDW=.TRUE.
AGGAC=0.0000
  • vdW-DF2:采用该范德瓦尔斯修正需在INCAR中加上以下开关:
{.line-numbers} GGA
1
2
3
LUSE_VDW=.TRUE.
Zab_vdW=-1.8867
AGGAC=0.0000

在考虑DFT-D2或者DFT-D3计算的过程中可同时考虑自旋轨道耦合(SOC),但是在vdW-DF、optB88-vdW 、optB86b-vdW和vdW-DF2中无法同时考虑自旋轨道耦合进行计算。

DFT+U

  1. 为什么加U?
    对于含有d、f轨道电子的强关联体系(Hubbard模型),电子间存在强烈的在位库仑相互作用,而交换相关泛函中的局域密度近似(LDA)或广义梯度近似(GGA)对电子之间的强在位库仑相互作用描述不准, 此时可将研究体系的交换相关泛函分成两部分计算,也就是用DFT+U的方法进行求解,公式如下:
    其中一部分电子用DFT算法(如LDA,GGA)等可以比较准确地描述;另外的d、f轨道电子通过引入Hubbard项得到正确的描述(U值就是考虑了同一个原子自旋相反的局域电子之间的库仑排斥)
  2. 如何加U?
    DFT计算时输入文件的U值由以下参数确定:
    {.line-numbers}
    1
    2
    3
    4
    5
    LDAU= .TRUE.|.FALSE.:开启/关闭+U功能,默认值为.FALSE.;
    LDAUTYPE=1|2|4指的是+U的类型,默认值是2;其中1Liechtenstein等提出的旋转不变LSDA+U方法;2Dudarev等提出的简化 LSDA+U方法;41类似, 但不考虑LSDA交换劈裂;
    LDAUL=-1|1|2|3分别对应不加U、p、d、f轨道加U;
    LDAUULDAUJ分别指定电子库伦相互作用项和交换相互作用项(U和J值);
    LMAXMIX =2/4/6:默认为2,加U计算时该值需大于轨道量子数,对于含有d轨道或f轨道电子的体系需对应增加至46
    注意:
    1)必须为每一类原子指定一个U值,其顺序与结构文件POSCAR中元素排序一致;
    2)同一种元素在不同的体系中U值不同;
    3)U在所有电子中都存在,不管自旋是否相同,J只存在于自旋相同的电子上;当LDAUTYPE =2时,U-J的差值才有意义,即有效的U参数。
    该参数需要测试或者参考文献。

磁性体系

很多体系都有磁性,磁性体系的计算需要打开自旋。以FeO为例:

{.line-numbers}
1
2
ISPIN=2
MAGMOM=原子数目*初始磁矩,eg 3*2 4*5 (注:*号前后不要有空格,且对应顺序)

对于一些简单的磁性体系,我们可以直接使用ISPIN=2, MAGMOM不必进行设置。对于复杂的体系则需要设置。初始值不需要与实验值完全一致,一般设置1.5倍即可。

  1. 有4个Fe原子,如果考虑铁磁态,则可以设置 MAGMON= 4*3 4*0 (Fe初设的磁矩为3,O元素为0)

  2. 有4个Fe原子,如果考虑反铁磁态,则可以设置 MAGMON=3 -3 3 -3 0 0 0 0或者MAGMON=3 3 -3 -3 0 0 0 0

表面优化

1. Bulk计算,获得稳定构型。

注意LREAL、PREC、ENCUT、EDIFF、ISMEAR保持一致。

2. 考虑不同晶面、层数。

对于表面结构,有以下几个需要注意的:

  • xy 方向上表面的大小(15 Å或者100原子左右);
    A 这个影响表面吸附物种的覆盖度;
    B 影响体系的尺寸大小和计算时间;
    C 不同的大小需要选取对应的K点;回想一下我们前面提到的经验规则。
  • 不同的晶面,(111), (100),(110);
    A 这取决于你研究的方向;低指数面、解离面。
    B 不同晶面的表面能;
    C 不同晶面的表面结构,反应活性等。
  • 表面结构的层数(3+3以上,6Å+6Å)
    A 层数多了,原子变多,体系在 z 方向的尺寸增加,也会影响计算速度;
    B 计算中需要弛豫的层数;
    C 不同层数对你要计算的性质会产生影响,比如表面能;
    D 不同晶面需要的层数也不尽相同,一般开放的表面需要更多的层数;
    E 根据自己感兴趣的性质,选择合适的层数,也就是需要你去测试一番。
  • Slab模型有两种,一种是上下表面对称的,一种是非对称的。对称性结构往往需要很多层,体系较大。 非对称的结构体系较小,但存在偶极矩的影响,要注意加LDIPOL 和IDIPOL这两个参数来消除。1、2、3分表代表在x、y和z 方向上进行校正,4代表在所有方向上。真空层15-20埃。

3. 工具:MS、ASE、P4vasp、vaspkit。

MS操作。

一般用conventional cell。
MS操作。Windows
VESTA将CONTCAR转成cif,导入MS。
build-surface-cleave surface 切面
Surface mesh-surface vectors 可调整vector。

Surface box-top/thickness,调整终端和厚度。
Build-crystals-build vaccum slab-图2-build。

Lattice parameters-advanced-re-orient to startard把真空层沿Z轴 图3。

转换成POSCAR,可优化后进一步扩胞。

ASE操作。python

块体优化后
module load python/3.7.3
python

from ase.io import read
from ase.io import write
from ase.build import surface
a=read(”CONTCAR”,format=’vasp’)
b=surface(a,(1,1,1),6,vacuum=7.5) 6层,真空层15Å
write(”POSCAR”,b,format=’vasp’)
Ctrl+D (退出 python)

大师兄ASE脚本:切金属

python3 cssm.py

vaspkit

块体优化后复制成POSCAR
vaspkit –task 803
生成SLAB***.vasp

4 表面弛豫计算

单点能

转换成POSCAR单点计算。或表面优化时第一个离子步的能量,但要保证这一步收敛。
KPOINTS设置,20-30。
INCAR设置,NSW=0或默认0.偶极矩矫正 LDIPOL = .TRUE. ; IDIPOL = 3。
POTCAR准备,提交任务脚本。完成后 grep without OUTCAR 获取单点能。

表面优化

1)POSCAR第7行后面加一行加S, Sed –i ‘7a S’ POSCAR
原子坐标加F/T。 Sed –i ’10,11s/$/ F F F/g’ POSCAR
Sed –i ’12,13,14s/$/ T T T/g’ POSCAR
vaspkit –task 402/403 或者ctrl+v,选中列-大写A-输入 T T T - ESC.
2)POTCAR和KPOINTS不变。
INCAR设置:IBRION=1、2,EDIFF=1E-6,EDIFFG=-0.03(比块体小一点),ISIF=2,
提交任务。

5 表面能计算

表面能是创造物质表面时,破坏分子间化学键所需消耗的能量。在固体物理理论中,表面原子比物质内部的原子具有更多的能量,因此,根据能量最低原理,原子会自发的趋于物质内部而不是表面。表面能的另一种定义是,材料表面相对于材料内部所多出的能量。把一个固体材料分解成小块需要破坏它内部的化学键,所以需要消耗能量。如果这个分解的过程是可逆的,那么把材料分解成小块所需要的能量就和小块材料表面所增加的能量相等。但事实上,只有在真空中刚刚形成的表面才符合上述能量守恒。因为新形成的表面是非常不稳定的,它们通过表面原子重组和相互间的反应,或者对周围其他分子或原子的吸附,从而使表面能量降低。
γ=($E_{surf}$ - N * $E_{bulkatom}$)/(2A)

1ev/$Å^2$ =16.02J/$m^2$
$E_{surf}$是未弛豫的表面能量,A是表面积a*b。

吸附能扫描

算吸附能时一般只需要算分子在高对称点的吸附能,比如算Li在铜表面的吸附能,我们只需算Li在Cu(111)面的Top,HCP,FCC和Bridge位的吸附能即可。但是这些离散的点没法帮你构建一个“势能面”,也就是无法得到吸附能在表面究竟是怎么分布的. 下面的文献里DOI:10.1063/1.4901055有一张图非常生动,讲的分别是Li在Li(001)以及Mg在Mg(0001)表面的吸附势能分布,我们不仅可以知道哪些位置有利于吸附,还可以根据吸附能分布,描绘出分子在表面的扩散路径,这有助于我们使用NEB方法算扩散能垒。

图片

其实思路很简单,就是将表面网格化,算分子在格点的吸附能,再画出等高线图。但是大部分吸附位点其实是不稳定的,所以我们采取固定吸附分子x,y方向,在z方向弛豫以达到吸附平衡的策略。

其实手动撒点,再采集数据也是可行的,但是会比较麻烦,因此作者根据实际需要开发了一款脚本scan_adsorption_energy用于自动完成这个过程。脚本使用Python编写,需要numpy和matplotlib第三方库。 我们首先算好一个吸附例子得到CONTCAR,这个可以让我们得到吸附分子的元素信息和理想的吸附高度。如下图,吸附的Li是第97号,也是最后一个原子,接着就可以变换这个结构得到不同的吸附结构。

图片

首先我们需要在目录下准备以下几个文件CONTCAR,INCAR,KPOINTS,POTCAR以及排队脚本文件(我的是vasp.pbs),然后使用脚本产生不同的结构:

tamas@tamas-PC:~/Desktop/scan$ python scan_adsorption_energy.py -g
Input which atom you want to move to scan adsorption energy–>-1
脚本提示我们移动第几个原子产生不同的吸附结构,我们可以输入97或者-1,移动最后一个Li。这个Cu(111)面是(331)的超胞,因此我们最需要取左下角的1/9的小区域撒点即可,然后程序会将小区域内的数据点进行平移得到整个表面的数据。如何根据不同的表面调整这个参数以及撒点的密度呢?可以使用-u参数进行自定义设置。

Input how many supercells in x axis–>3
Input how many supercells in y axis–>3
Input interporate how many points for each unit cell in x axis,even number are suggested–>2
Input interporate how many points for each unit cell in y axis,even number are suggested–>2
command to submit jobs,e.g. qsub vasp.pbs–>qsub vasp.pbs
Input which atom you want to move to scan adsorption energy–>-1
手动选择x方向的超胞数3,y方向的超胞数3,在小区域内的x方向的撒点数,y方向的撒点数(为了包含中心点,建议选择偶数)和提交排队脚本的命令(我的是qsub vasp.pbs)。脚本会自动产生4个结构(2*2),并会将每个计算所需的文件都复制到相应的文件夹内,并提交排队脚本进行计算。同时脚本会写出一个grid.info里面包含了超胞数,选择移动的原子以及晶格常数信息。 在所有的计算完成后,同样使用该脚本提取数据。

python scan_adsorption_energy.py -e
脚本会自动从每个结构文件夹内读取能量数据和吸附原子的坐标。生成两个dat文件new_direct.dat和new_dat-3D.dat。前者是原子的分数坐标和总能量,后者有4列,分别为原子的笛卡尔x分量坐标,y分量坐标,z分量坐标以及总能量。我们需要使用第1,2,4列绘制等高线图。如果需要绘制吸附能,需要自己手动将第四列的数据减去吸附原子和表面的能量。 脚本会尝试使用matplotlib绘制等高线图,如下图,你也可以使用Origin Pro绘图。

图片

脚本可以在我的Github仓库下载(https://github.com/tamaswells/VASP_script/blob/master/scan_adsorption_energy.py).

固定基矢优化结构

在第一性原理计算时,二维材料的建模会采用在某一方向(一般是z轴方向)添加15-20 Å的真空层的方式来抵消层与层之间的相互作用。但是,当我们优化晶格常数的时候会遇到一个问题:那就是使用ISIF=3来优化平面内晶格常数(一般是a和b)时,c轴方向的晶格长度总是随着优化变的越来越小。

  1. 采用ISIF=4的固定晶胞体积的方法来计算。这样采用实验参数建的晶胞,由于a和b方向的变化,c轴也会变化。解决了越是优化c轴越小的问题,但是计算量还是很大。

  2. 编译vasp,就自己编译了固定z轴的脚本,专门计算二维材料。这种方法优点是直接换脚本就能解决上面优化时c轴缩小的问题,但是我想固定两个轴就需要再编译一份脚本。比较麻烦。

修改constr_cell_relax.F实现固定基矢优化结构的原理:

FCELL contains the forces on the basis vectors. These forces are used to modify the basis vectors according to the following equations:

FCELL 是个3 * 3的矩阵,包含了晶格矢量的受力,当ISIF = 3时,晶格矢量会按照如下的方法更新:

A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style
DO J=1,3
DO I=1,3
DO K=1,3
A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
ENDDO
ENDDO
ENDDO
当通过OPTCELL文件把晶格矢量在某个方向的受力改成0以后,

IF (ICELL(I,J)==0) FCELL(I,J)=0.0

在结构更新的时候,该方向就不会改变。完成编译后,在 INCAR 文件中设置 ISIF=3,便可实现固定 Z 轴的结构优化。
对于其他方向晶格基矢的修改同理:对于a方向基矢,将 FCELL(3,I) 和 FCELL(I,3) 分别改为 FCELL(1,I) 和 FCELL(I,1) ;对于b方向基矢,则分别改为 FCELL(2,I) 和 FCELL(I,2) ;固定两个基矢则应该同时将两个方向对应的矩阵元设置为0。或参考 VASP并行可执行软件包,可对晶胞参数进行部分优化 可实现添加参数文件实现不同基矢方向固定。(未做测试)

  1. 肖承诚博士的补丁,按照肖博士的方法安装以后,在INCAR中添加如下参数即可。

IOPTCELL = xx xy xz yx yy yz zx zy zz

上面的九个参数对应晶格常数矩阵。

1 xx xy xz

2 yx yy yz

3 zx zy zz

1代表驰豫,0代表固定

此时固定c轴的设置就是如下:

IOPTCELL = 1 1 0 1 1 0 0 0 0

下载地址:https://github.com/Chengcheng-Xiao/VASP_OPT_AXIS

滑移能和剥离能

材料的塑性行为主要是由于材料的结构和化学键决定的。材料发生塑性形
变在原子尺度上可以分解成两个过程。

一个是由于原子、位错或者晶界沿着某
一特定方向自由运动;第二个是在这种运动过程中存在束缚力保证在其它方向
不会产生裂纹,导致材料的断裂。因此任何一个塑性材料必须满足两个要求,
一是对某些原子层存在势垒相对低的滑动方式,以确保原子、位错或者晶界可
以沿着这些方向滑动。另一个是在这些滑动过程中,滑移平面两侧的原子依然
存在某种相对强的作用力,以确保滑移过程中材料的完整性。如金属键,其自
由电子容易发生变形和转移,在滑移过程中金属键不断的断裂和产生,从而导致了金属的滑移能量很小。而金属键的键强度很大,保证的金属滑移过程的完
整性。而离子键和共价键由于配位数和方向性等要求,在滑移时必须先发生键
的断裂,从而阻碍滑移的进行。由范德华力连接的层状材料,虽然层与层之间
容易发生滑移,但其层间作用力较弱,材料在滑移过程中容易发生破坏。

为了研究 α-Ag2S 的塑性行为,我们通过第一性原理对 α-Ag2S 的滑移过程进
行了计算。其中,我们计算的滑移面为(100)方向,滑移方向为[001]方向。其
滑移过程如图 6.10 所示,α-Ag2S 的滑移面沿着能量最低点进行滑移,滑移面不
仅发生相对移动,其面间距也会发生。因此,整个滑移路径成波浪形。 我们将滑移一个周期(滑移一个晶格常数)均分 12 步。每一步我们计算系
统总能量与滑移面之间距离的关系,对于塑性材料而
言,它需要同时具有足够小的滑移能来保证材料的滑移可以进行;以及大的解离
能以至于材料不会发生断裂。

这里主要分享怎样快速开始第一个结构优化过程。

Linux常用命令

目录操作

1 切换目录(cd)

{.line-numbers}
1
2
3
4
5
6
7
cd /                 //切换到根目录
cd /bin //切换到根目录下的bin目录
cd ../ //切换到上一级目录 或者使用命令:cd ..
cd ~ //切换到home目录
cd - //切换到上次访问的目录
cd xx(文件夹名) //切换到本目录下的名为xx的文件目录,如果目录不存在报错
cd /xxx/xx/x //可以输入完整的路径,直接切换到目标目录,输入过程中可以使用tab键快速补全

2 查看目录(ls)

{.line-numbers}
1
2
3
4
ls                   //查看当前目录下的所有目录和文件
ls -a //查看当前目录下的所有目录和文件(包括隐藏的文件)
ls -l //列表查看当前目录下的所有目录和文件(列表查看,显示更多信息),与命令"ll"效果一样
ls /bin //查看指定目录下的所有目录和文件

3 创建目录(mkdir)

{.line-numbers}
1
2
mkdir tools          //在当前目录下创建一个名为tools的目录
mkdir /bin/tools //在指定目录下创建一个名为tools的目录

4 删除目录与文件(rm)

{.line-numbers}
1
2
3
4
5
6
rm 文件名              //删除当前目录下的文件
rm -f 文件名 //删除当前目录的的文件(不询问)
rm -r 文件夹名 //递归删除当前目录下此名的目录
rm -rf 文件夹名 //递归删除当前目录下此名的目录(不询问)
rm -rf * //将当前目录下的所有目录和文件全部删除
rm -rf /* //将根目录下的所有文件全部删除【慎用!相当于格式化系统】

5 修改目录(mv)

{.line-numbers}
1
2
3
mv 当前目录名 新目录名        //修改目录名,同样适用与文件操作
mv /usr/tmp/tool /opt //将/usr/tmp目录下的tool目录剪切到 /opt目录下面
mv -r /usr/tmp/tool /opt //递归剪切目录中所有文件和文件夹

6 拷贝目录(cp)

{.line-numbers}
1
2
cp /usr/tmp/tool /opt       //将/usr/tmp目录下的tool目录复制到 /opt目录下面
cp -r /usr/tmp/tool /opt //递归剪复制目录中所有文件和文件夹

7 搜索目录(find)

{.line-numbers}
1
find /bin -name 'a*'        //查找/bin目录下的所有以a开头的文件或者目录

8 查看当前目录(pwd)

{.line-numbers}
1
pwd                         //显示当前位置路径

文件操作

1 新增文件(touch)

{.line-numbers}
1
touch  a.txt         //在当前目录下创建名为a的txt文件(文件不存在),如果文件存在,将文件时间属性修改为当前系统时间

2 删除文件(rm)

{.line-numbers}rm 文件名link
1
rm -f 文件名           //删除当前目录的的文件(不询问)

3 编辑文件(vi、vim)

{.line-numbers}
1
vi 文件名              //打开需要编辑的文件

–进入后,操作界面有三种模式:命令模式(command mode)、插入模式(Insert mode)和底行模式(last line mode)
命令模式
-刚进入文件就是命令模式,通过方向键控制光标位置,
-使用命令”dd”删除当前整行
-使用命令”/字段”进行查找
-按”i”在光标所在字符前开始插入
-按”a”在光标所在字符后开始插入
-按”o”在光标所在行的下面另起一新行插入
-按”:”进入底行模式
插入模式
-此时可以对文件内容进行编辑,左下角会显示 “– 插入 –””
-按”ESC”进入底行模式
底行模式
-退出编辑: :q
-强制退出: :q!
-保存并退出: :wq

操作步骤示例

1.保存文件:按”ESC” -> 输入”:” -> 输入”wq”,回车 //保存并退出编辑
2.取消操作:按”ESC” -> 输入”:” -> 输入”q!”,回车 //撤销本次修改并退出编辑

补充

{.line-numbers}
1
2
vim +10 filename.txt                   //打开文件并跳到第10行
vim -R /etc/passwd //以只读模式打开文件

4 查看文件

{.line-numbers}
1
2
3
4
cat a.txt          //查看文件最后一屏内容
less a.txt //PgUp向上翻页,PgDn向下翻页,"q"退出查看
more a.txt //显示百分比,回车查看下一行,空格查看下一页,"q"退出查看
tail -100 a.txt //查看文件的后100行,"Ctrl+C"退出查看

5 grep

{.line-numbers}
1
2
3
grep -i "the" demo_file              //在文件中查找字符串(不区分大小写)
grep -A 3 -i "example" demo_text //输出成功匹配的行,以及该行之后的三行
grep -r "ramesh" * //在一个文件夹中递归查询包含指定字符串的文

文件权限

1 权限说明

文件权限简介:’r’ 代表可读(4),’w’ 代表可写(2),’x’ 代表执行权限(1),括号内代表”8421法”
##文件权限信息示例:-rwxrw-r–
-第一位:’-‘就代表是文件,’d’代表是文件夹
-第一组三位:拥有者的权限
-第二组三位:拥有者所在的组,组员的权限
-第三组三位:代表的是其他用户的权限

2 文件权限

普通授权 chmod +x a.txt
8421法 chmod 777 a.txt //1+2+4=7,"7"说明授予所有权限

打包与解压

1 说明

{.line-numbers}
1
2
3
4
.zip、.rar        //windows系统中压缩文件的扩展名
.tar //Linux中打包文件的扩展名
.gz //Linux中压缩文件的扩展名
.tar.gz //Linux中打包并压缩文件的扩展名

2 打包文件

tar -zcvf 打包压缩后的文件名 要打包的文件
参数说明:z:调用gzip压缩命令进行压缩; c:打包文件; v:显示运行过程; f:指定文件名;
示例:

{.line-numbers}
1
tar -zcvf a.tar file1 file2,...      //多个文件压缩打包

3 解压文件

{.line-numbers}
1
2
3
4
tar -zxvf a.tar                      //解包至当前目录
tar -zxvf a.tar -C /usr------ //指定解压的位置
unzip test.zip //解压*.zip文件
unzip -l test.zip //查看*.zip文件的内容

原文链接:https://blog.csdn.net/m0_46422300/article/details/104645072

其他准备

安装好VASP、vaspkit,个人电脑安装好VESTA、MS或者qvasp等建模可视化软件中的一个。如果是租用超算,还需要准备好提交任务的脚本文件,一般工程师会帮我们编好,软件安装只需要给他们提供安装包或官网即可。

vaspkit运行的五种方式

方法一 :交换式菜单模式,直接输入vaspkit;

方法二 :命令+参数模式,适合批处理。如vaspkit -task 102 -kpr 0.03,详细介绍见vaspki-help, 注意部分功能不支持这种模式;

方法三 :利用bash管道传递参数,如echo -e “102\n2\n0.04\n”| vaspkit

方法四 :同样利用bash管道传递参数, (echo 102; echo 2; echo 0.04)|vaspkit

方法五 :vi cmd.in(文件名任意,非必须是cmd.in)包含以下三行内容

{.line-numbers}
1
2
3
102
2
0.04

然后运行vaspkit < cmd.in,产生KPOINTS文件,其他命令格式类似。

输入文件

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
INCAR      in    **
STOPCAR in
stout out
POTCAR in **
KPOINTS in **
IBZKPT out
POSCAR in **
CONTCAR out
CHGCAR in/out
CHG out
WAVECAR in/out
TMPCAR in/out
EIGENVAL out
DOSCAR out
PROCAR out
OSZICAR out
PCDAT out
XDATCAR out
LOCPOT out
ELFCAR out
PROOUT out

  • 必要的输入文件有四个:INCAR,KPOINTS,POSCAR,POTCAR
    • INCAR 告诉VASP算什么,怎么算。
    • KPOINTS 包含计算的倒空间K点信息。
    • POSCAR 是结构信息,包括元素、晶胞参数及各个原子的坐标信息。
    • POTCAR 元素的赝势文件,描述体系中对应的原子核和电子的相关信息。
    • 提交任务的脚本或者命令,需要你自己准备。

下面由简单到复杂介绍

POSCAR

结构文件信息格式如下:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
Mo S2                  #结构名称,可不写或随便写,不能没有这一行
1.0 # universal scaling parameters,缩放系数
3.1500000954 0.0000000000 0.0000000000 # 晶格矢量 a(1)
-1.5750000477 2.7279801045 0.0000000000 # 晶格矢量 a(2)
0.0000000000 0.0000000000 12.3000001907 # 晶格矢量 a(3)
Mo S # 元素名称,种类
2 4 #对应元素的原子个数,顺序要对应
Direct #
0.333333345 0.666666689 0.250000000
0.666666624 0.333333323 0.750000019
0.333333345 0.666666689 0.620999986
0.666666624 0.333333323 0.379000014
0.666666624 0.333333323 0.121000053
0.333333345 0.666666689 0.878999976
  • 对于自己建模的结构,按照该格式自己写;
  • 对于Materials projiect数据库可以直接下载POSCAR格式;
  • 对于.cif文件,可以导入 VEASTA 等输出为POSCAR的格式,或者用 vaspkit 的411/412功能。

Mind:

POSCAR文件非常重要,大部分报错都是结构不合理导致的。一个好的初始结构不仅可以加快收敛速度,也会增加准确性。文件中元素的顺序非常重要。

KPOINTS

文件一般格式如下:

{.line-numbers}
1
2
3
4
5
K-POINTS      #  第一行随便写或者空着不写,但不能没有这一行
0 # 零,表示自动生成
Gamma # gamma点centered
1 1 1 # 1*1*1格子,K点密度
0 0 0 # S1 S2 S3,shift的值,一般保持 0 0 0 不变。

大家需要修改的是第三行和第四行,其他可不改。

  • 第三行:VASP只认第一个字母,大小写均可。当然这一行也可以直接写字母G或者g。G表示的是以gamma点为中心生成网格。
    另外一种是原始的Monkhorst-Pack 网格,两者的区别是 M 或者 m 在 G的基础上在三个方向上平移了1/(2N)个单位。G,也叫 gamma centered Monkhorst-Pack Grid;所以,gamma centered 只是MP网格的一种特殊情况。
    因此这一行也可以一直用G不变。

  • 第四行,在xyz三个方向上生成对应数目的K点;__这个一般需要测试__。如果开始不想测试,官网也给出了经验方法,

    • d区金属,k*a ~ 30 Å
    • 普通金属,k*a ~ 25 Å
    • 半导体,k*a ~ 20 Å
    • 绝缘体,k*a ~ 15 Å

    其中a对应三个方向的晶胞参数,也就是对应方向的k点与该方向晶胞参数乘积最好在20-30之间。

如果已经安装了vaspkit程序,可以通过102功能自动生成KPOINTS文件。

POTCAR

POTCAR赝势文件相对来说最复杂也最简单。其中有很多信息,大部分信息我们不需要关注,一般只需要关注前几行。特别是赝势的种类信息。以Mo为例:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
PAW_PBE Mo 08Apr2002                   
6.00000000000000
parameters from PSCTR are:
VRHFIN =Mo: 4p5s4d
LEXCH = PE
EATOM = 217.5176 eV, 15.9871 Ry

TITEL = PAW_PBE Mo 08Apr2002
LULTRA = F use ultrasoft PP ?
IUNSCR = 1 unscreen: 0-lin 1-nonlin 2-no
RPACOR = 2.200 partial core radius
POMASS = 95.940; ZVAL = 6.000 mass and valenz
RCORE = 2.750 outmost cutoff radius
RWIGS = 2.750; RWIGS = 1.455 wigner-seitz radius (au A)
ENMAX = 224.584; ENMIN = 168.438 eV

除此之外,一种元素也有很多赝势,例如Mo:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
[scf1355@ln151%nc-l potpaw_PBE.52]$ ls Mo*
Mo:
POTCAR PSCTR

Mo_pv:
POTCAR PSCTR

Mo_sv:
POTCAR PSCTR

Mo_sv_GW:
POTCAR PSCTR

如果没有特别的需求,直接采用VASP官网推荐的即可。参考链接:
https://cms.mpi.univie.ac.at/vasp/vasp/Recommended_PAW_potentials_DFT_calculations_using_vasp_5_2.html
不知道选哪一个就选第一个,使用方法如下:

{.line-numbers}
1
cat ~/pot/Mo/POTCAR ~/pot/S/POTCAR  > POTCAR

或者用vaspkit的103命令自动生成。

Mind:

该文件看似复杂,使用时只需要查看元素种类和POSCAR是否对应即可。
可以通过grep PBE POTCAR查看顺序是否正确。

{.line-numbers}
1
2
3
4
5
[scf1355@ln151%nc-l 2h]$ grep PBE POTCAR 
PAW_PBE Mo 08Apr2002
TITEL = PAW_PBE Mo 08Apr2002
PAW_PBE S 06Sep2000
TITEL = PAW_PBE S 06Sep2000

INCAR

这个控制计算过程的文件最为复杂,手册上相关参数多达三百多个。因此前期我们记住一点:不知道含义的参数不要加。不加对于有些参数意味着用默认值,一般是合理的。INCAR参数查找参考官方:https://www.vasp.at/wiki/index.php/Category:INCAR_tag

贴一个我的结构优化的INCAR:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ISTART =  0            (Read existing wavefunction; if there)
ISPIN = 2 (Non-Spin polarised DFT)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = A (Precision level)
LWAVE = .F. (Write WAVECAR or not)
LCHARG = .F. (Write CHGCAR or not)
ISMEAR = 0 (Gaussian smearing; metals:1)
SIGMA = 0.05 (Smearing value in eV; metals:0.2
EDIFF = 1E-08 (SCF energy convergence; in eV)
EDIFFG = -1E-02 (Ionic convergence; eV/AA)
NSW = 100 (Max ionic steps)
IBRION = 2 (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
POTIM = 0.1
ISIF = 3 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)

上面参数详细介绍

  • ISTART :初始波函数。WAVECAR存在默认1,不存在默认0。
    0:随机生成波函数
    1:从 WAVECAR 读取波函数,若读取失败则随机生成。
  • ICHARG :初始电荷密度。当ISTART=0默认2,其他默认1。
    0:通过 WAVECAR 生成电荷密度,若读取失败则通过原子电荷密度叠加生成
    1:从 CHGCAR 读取电荷密度
    2:通过原子电荷密度叠加生成电荷密度
    11:从 CHGCAR 读取电荷密度,且在(非自洽)循环中保持不变
  • ISPIN :自旋极化。默认1,关闭。不知道是否打开,可以打开,结果不会有影响。
    1:非自旋极化(每个轨道上自旋向上和自旋向下的电子数量相等,适合非磁性体系)
    2:自旋极化(适合铁磁、反铁磁体系)
  • ENCUT :平面波截断能(单位 eV)
    默认元素中最大的ENMAX。一般为 POTCAR 中 最大ENMAX 参数值的 1.0 到 1.3 倍。500适用大部分体系。
  • PREC :总体计算精度,默认Normal。可设为 Low、Med、High、Normal、Single、Accurate。具体参照手册,设为A或者Accurate足够精确。
  • LWAVE LCHARG :控制WAVECAR和CHGCAR输出,默认输出,在结构优化可以关闭,节省储存空间。
  • ISMEAR :轨道分数占据
    −5:四面体方法(适合半导体和绝缘体)
    0:高斯方法(适合导体、半导体和绝缘体),展宽由 SIGMA 确定
    N:Methfessel-Paxton 方法(适合导体)
  • EDIFF :自洽循环收敛标准(系统能量变化),单位为 eV。默认1E-4,可以设的更精确一些。一般设置为 1E-6或更高。
  • __EDIFFG__:离子位置优化收敛标准。默认EDIFF*10。正值为系统能量变化,负值为原子上残余力。随着体系维数下降,可适当增加原子残余力。
  • NSW :离子位置优化最大步数(IBRION=1、2)。默认0。设为100足够大,也可以继续加大。
    分子动力学开启时代表分子动力学模拟步数(IBRION=0)
  • IBRION :离子位置优化算法。默认NSW=0/1时为-1,其他默认0。
    -1:离子不移动
    0:分子动力学模拟
    1:准牛顿法
    2:共轭梯度法
    5:振动频率计算
    6:弹性常数计算
    7/8:DFPT方法
  • POTIM :有限差分法步长。当IBRION=1/2/3默认值 0.5。当IBRION=0(MD)时必须指定。
  • ISIF :需要结构优化变量。默认IBRION=0时为0,其他默认2。结构优化一般设为3,全部优化。

Mind

  • 上面四个输入文件都是纯文本文件,且文件名不能更改,格式不能错误。
  • vaspkit可以同时产生INCAR、POTCAR、KPOINTS文件。

提交任务

在超算上提交或者集群上提交,只需要把四个输入文件和提交任务的脚本文件放在同一文件夹下,运行脚本即可。
在个人服务器、虚拟机提交可以通过nohup mpirun -np 核数 vasp_std > std.out &提交

Mind

  • 提交后经常查看作业收敛情况,有问题及时暂停,避免浪费时间

查看任务完成

实时查看任务日志tail -f 日志文件
任务停止后,判断任务是否正常结束,输入grep reached 日志文件或者tail 日志文件
如果输出reached required accuracy - stopping structural energy minimisation代表正常结束,否则很可能没有收敛或者出错中断。

输出文件

任务一旦提交,会立刻生成输出文件,不过此时可能是空的。运行片刻大部分文件会产生内容。
https://www.vasp.at/wiki/index.php/Category:Output_files
其中我们先关注几个输出文件。

CONTCAR

优化后的结构文件,与POSCAR相对应。为最后一个离子步的结构。当然如果要提高精度或者续算,可以cp CONTCAR POSCAR继续计算。

OSZICAR

包含自洽计算中能量收敛等信息。

OUTCAR

这个是VASP最主要的输出文件,包含计算过程中大量信息,依次主要包括:
VASP版本;
计算开始时间和并行性CPU数;
赝势信息;
最近邻列表;
对称性信息;
晶格信息和k点坐标;
INCAR中读入的参数和其他大部分的默认参数值;
平面波个数和FFT信息;
每一步离子步数和其中每一个电子自洽的时间、内存、能量等信息;
自洽完成后的费米能和能量本征值;
应力;力;
电荷数和磁矩;
程序运行时间。
总而言之想找而不知道去哪里找,就在这里。

IBZKPT

由VASP自动生成,包含所有不等价k点的坐标和权重以及可能的四面体链接情况,它的格式与KPOINTS文件的格式是完全一样的。

WAVECAR

波函数文件,二进制文件,不能直接用文本编辑器直接打开。
波函数文件较大,可以通过(LWAVE)来控制输出。结构优化不需要输出。
它也可以作为输入文件,为后续计算提供初始波函数。

CHGCAR和CHG

电荷密度文件,包含晶格矢量,原子坐标,总电荷密度以及PAW的单中心占据情况等信息。 该文件既是输出文件(由LCHARG控制),也可以是输入文件。结构优化不需要输出。
在计算能带或者态密度时,需要读入电荷密度进行非自洽计算。

总结

  • 输入文件KPOINTS比较简单,可以手动,注意格式;INCAR参数众多,直接参考,不知道的参数不加;POTCAR文件比较大,我们只关注赝势对应的元素顺序与POSCAR是否一致;POSCAR结构文件VESTA直接输出。提交任务脚本由超算代理提供,不同超算不同。
  • 任务提交后关注是否正常结束,另外关注计算的时间等。
0%