2024年2月22日发(作者:)
使用VASP的个人经验手册
侯柱锋
厦门大学物理系2004届博士
E-mail: zhufhou@
2004/06/22
本手册纯属个人使用VASP后的心得和经验总结,版权属于本手册的作者及厦门大学物理系计算物理实验室(Group leader: 朱梓忠教授)。未经许可,不准在网上传阅。文中提到的一些小程序,可以提供使用。在参考的过程,如遇到不清楚或含糊的地方,可以参考VASP的英文manual或email给我。如认为本手册某些地方需要更正或修改的,请email给我。当在使用VASP的过程中遇到问题,也可以email给我,大家一起学习VASP的使用,挖掘和掌握VASP强大的功能。本手册参考了VASP的英文manual、的报告以及从internet网上收集的资料。
本手册大致有以下几个内容:
A程序的编译
“ñAVASP的主要输入文件
OAVASP的主要输出文件
lA参数设置与选择的技巧
A材料基态性质的计算方法和步骤
ZA材料磁性性质的计算
µA表面体系的计算
”ªAtools中小程序的说明
A半导体中的缺陷和杂质问题(暂未完成)
十、如何进行分子动力学模拟(暂未完成)
十一、强关联体系的计算(LDA+U或GGA+U)(暂未完成)
一、程序的编译
声明:本实验室(厦门大学物理系计算物理实验室, Group leader: 朱梓忠教授)购买的是VASP4.4.5版本,所属本实验室的成员以及经过朱梓忠教授同意使用的合作者必须遵守该软件的使用协议,注意VASP软件的版权问题,严禁私下发布或传播本实验室购买的VASP源代码和赝势库以及编译VASP得到的可执行代码。
1
下面以编译VASP4.4.5版本为例,编译更新的版本VASP4.5.5、VASP4.6和VASP5.0(即将发布)的步骤与此相同。
1、 所需文件和程序
VASP源代码:和
数学库:LAPACK和BLAS (/),
或mkl(配合intel的fotran编译器用),
或ATLAS (/)
或Lib GOTO (/users/flame/goto/)
Fortran编译器:PGI fortran 至少4.0以上版本(/),
或Intel的 ifc (8.0以上版本是ifort, /software/products/compilers/flin/),前者可以从网站上下载到15天的试用版本,后者可以从网站下载到免费的版本。或者在国内的个人ftp服务器上搜索它们的破解版本。
本实验室的都有这些软件的备份。
2、下面采用PGI fortan编译器pgf90、ATLAS数学库对VASP4.4.5进行编译
这里假定已经安装好了fortran编译器,所有文件都放在/home/houzf/VASP_SRC目录下,机器的操作系统是Linux: Redhat9.0。
a) 从/下载atlas3.6.0_Linux_,并用如下命令解压:tar xzvf atlas3.6.0_Linux_
解压后得到一个目录Linux_P4SSE2,在此目录下有个lib子目录,该lib子目录中的文件为libatlas.a, libcblas.a, libf77blas.a, liblapack.a, 这些就是编译vasp时所需要的数学库文件之一。
b) 用如下命令解压和:
tar xzvf
tar xzvf
解压后分别得到目录vasp.4.4和,目录vasp.4.4中文件是vasp的主要源代码,是编译vasp时需要的一些特定的数学库程序,在这两个目录中都有编译时所用的makefile文件,针对机器和fortran编译器,选择相应的makefile。
c) 进入目录,选择_pg,并把它拷贝成makefile,然后键入make命令开始编译。整个命令如下:
cd
cp _pg makefile
make
编译成功后,得到libdmy.a文件。
d) 退出目录,进入vasp.4.4目录,选择_pg,并把它拷贝成makefile,编辑makefile文件,通过修改LIB变量的赋值而采用基于ATLAS的数学库文件,修改的地
2
方和方法是:
在第87和88行前加上#,把这两行注释掉,然后去掉第91,92和93行前的#。
修改前和后的内容为分别为:
LIB = -L../ -ldmy ..//linpack_double.o
..//lapack_double.o -L/usr/local/lib /usr/local/lib/libblas.a
#
# the following lines should allow you to link to atlas based blas
#LIB = -L../ -ldmy ..//linpack_double.o
# ..//lapack_double.o -L/usr/local/lib
# -L$(HOME)/archives/BLAS_OPT/ATLAS/lib/Linux_ATHLONTB/ -lf77blas –latlas
#LIB = -L../ -ldmy ..//linpack_double.o
# ..//lapack_double.o -L/usr/local/lib /usr/local/lib/libblas.a
#
# the following lines should allow you to link to atlas based blas
LIB = -L../ -ldmy ..//linpack_double.o
..//lapack_double.o -L/usr/local/lib
-L../Linux_P4SSE2/lib/ -lf77blas -latlas
修改后保存makefile文件,键入make命令开始编译vasp。整个命令为:
cd ..
cd vasp.4.4
cp _pg makefile
编辑修改makefile文件
make
编译成功后,就可以得到VASP的可执行文件vasp。
e) 以root帐号登录机器,把成功编译VASP后得到的vasp放到/bin目录下,则任何一个普通用户都可以使用vasp。此时vasp可以当成于一个linux的命令来使用了,不再需要把vasp拷贝到当前的计算目录下。
二、VASP的主要输入文件
VASP的主要输入文件有INCAR, POTCAR, POSCAR和KPOINTS。INCAR文件控制了vasp进行何种性质的计算,POTCAR文件包含了体系中各类元素的赝势,POSCAR文件描述了所计算的体系的晶胞参数(包括基矢或平移矢量,晶格常数,原子位置等信息),KPOINTS描述了不可约布里渊区中k点取样,即k点设置。
1、INCAR文件
此文件控制vasp进行何种性质的计算,以及设置了计算方法中一些重要的参数。其中的关
3
键词可以分为如下几类:
对所计算的体系进行注释:SYSTEM
定义如何输入或构造初始的电荷密度和波函数:ISTART, ICHARG, INIWAV
定义价电子部分的如何驰豫:
平面波切断动能和缀加电荷时的切断值:ENCUT, ENAUG
电子部分优化的方法:ALGO, IALGO, LDIAG
电荷密度混合的方法:IMIX, AMIX, AMIN, BMIX, AMIX_MAG, BMIX_MAG, WC,
INIMIX, MIXPRE, MAXMIX
自洽迭代步数和收敛标准:NELM, NELMIN, NELMDL, EDIFF
定义离子芯部分的如何驰豫:
离子如何移动以及步长和步数:IBRION, NFREE, POTIM, NSW
分子动力学相关参数:SMASS, TEBEG, TEEND, POMASS,NBLOCK, KBLOCK,
PSTRESS
离子驰豫收敛标准:EDIFFG
定义态密度积分的方法和参数:
smearing方法和参数:ISMEAR, SIGMA
计算态密度时能量范围和点数:EMIN, EMAX, NEDOS
计算分波态密度的参数:RWIGS, LORBIT
其他:
计算精度控制:PREC
磁性计算:ISPIN, MAGMOM, NUPDOWN
交换关联函数:GGA, VOSKOWN
计算ELF和总的局域势:LELF, LVTOT
结构优化参数:ISIF
一般要设置的关键词:SYSTEM, ENCUT, ISTART, ICHARG, PREC, ISMEAR, SIGMA。针对计算不同的性质,再另外增加相应的关键词。
例子:
General:
SYSTEM = fcc Si !自洽计算fcc结构的Si
ISTART = 0 !开始新的计算
ICHARG = 2 !从原子的电荷密度重叠构造初始电荷密度
ENCUT = 240 !平面波切断动能
ISMEAR = 0; SIGMA = 0.1 !采用Gaussian smearing方法,展宽为0.1eV
PREC = Accurate !计算精度
4
2. POTCAR文件
赝势文件,最重要的输入文件之一。赝势库中赝势文件可以进行如下分类:
根据方法不同有Ultra-soft赝势(USPP)和投影缀加波的赝势(PAW)
根据交换关联函数的不同有LDA和GGA(又可以再分为PW91和PBE)
根据处理了半芯态有A, A_sv和A_pv的不同
根据ENMAX的大小有A, A_s和A_h的不同
如何准备?
如果你拿到的赝势文件的格式用相应的命令把各元素的赝势合并到一个文件POTCAR中:
a) 是以Z为扩展名的文件,用命令: zcat POTCAR.Z >>aa
b) 是解压后的文件POTCAR,用命令:cat POTCAR >>aa
(当有多类原子时,按POSCAR文件各类原子的顺序,依次使用上面的命令,把相应原子的POTCAR.Z合并到aa文件中)
c) 然后把aa文件移到到要计算的目录中(mv aa 计算的目录/POTCAR).
注释:在处理磁性材料,所计算的体系含有碱金属、碱土金属、周期表左边的3d过渡元素、镧系和锕系元素时,强烈推荐用PAW势,计算精度有提高。在采用超越赝势(USPP)时,使用PW91的GGA时,强烈要求把VOSKOWN = 1给选上。在采用PAW势时,一般推荐用LDA和PBE的。
下面给出PAW对不同元素,采用何种类型的PAW以及ENCUT值至少要取多少,所列的表格,供选择赝势时作为参考(下面几个表格中,红色表示是一般情况下首选用这种类型的PAW势,表格中数字表示的是切断动能值):
B_h 700
B 318
B_s 250
Al 240
Al_h 295
Ga 134
Ga_d 282
Ga_h 404
In 95
In_d 239
Tl 95
C_h 700 N_h 700O_h 700C 400
C_s 273
N 400N_s 250O 400F_h 700
F 400
O_s 250 F_s 250
S 280 Cl 280
S_h 402 Cl_h 409
Se 211 Br 216
Si 245 P 270Si_h 380 P_h 390Ge 173
Ge_d 287
Ge_h 410
Sn 103Sn_d 241Pb 98Sb 172
As 208Te 174 I 175
Bi 105
5
Tl_d 239
Pb_d 237Bi_d 242
注释:X_d表示的是,d电子作为半芯态来处理的。为了得到较高的计算精度,一般推荐采用X_d的赝势。X_h表示该势比较硬,也是切断动能要用的很大,它们一般是用含有这类原子的氧化物的计算中,为了提高计算的精度。其中Si_h一般用在含Si的沸石材料中。
H 250
H_h 700
Li 140 Be 300
Li_sv 271 Be_sv 308
Na 81
Na_pv 300
Na_sv 700
K_pv 150 Ca_pv 150K_sv 259 Ca_sv 290
Rb_pv 121 Sr_sv 226
Rb_sv 220
Cs_sv 220 Ba_sv 187
注释:这些元素一般很难赝化的,特别是与电负性很强的元素(比如F)结合时,计算的误差都比较大。X_sv表示把s电子作为半芯态处理,X_pv考虑把p电子作为半芯态来处理。
Mg 210
Mg_pv 265
Ti 178
Hf 220
Hf_pv 220
V 192
V_pv 263
Ta 223
Ta_pv 223
Ni 269
Ni_pv 367
Cr 227
Cr_pv 265
Mo 224
W 223
W_pv 223
Cu 273
Cu_pv 368
Ag 249
Au 229
Mn 269
Mn_pv 269
Tc 228
Tc_pv 228
Re 226
Re_pv 226
Zn 276
Cd 274
Hg
233
Sc_sv 222 Ti_pv 222
Y_sv 211 Zr_sv 229 Nb_pv 207 Mo_pv 224
Fe 267 Co 267
Fe_pv 293
Ru_pv 230
Os_pv 228
Ru 213 Rh 228 Pd 250
Rh_pv 271 Pd_pv 350
Pt 230 Os 228 Ir 210
6
注释:选择使用X_pv、X_sv还是X的赝势,一个与你要得到的计算精度有关,另外对这些元素在选择要注意些:
3d元素,一般选用X_pv,但是X的赝势也是能给出比较合理的结果。
4d元素,是最有问题的,强烈推荐用X_sv和X_pv的赝势。
5d元素,由于5p电子局域化很强,从Hf元素开始,可以选用X的赝势,推荐选用不同的赝势,进行test一下,然后选用合适的赝势。
Ce 300 Pr 252
La 219 Ac 169 Th 247 Pa 252
La_s 136
注释:如果f电子是itinerant(巡回的),则可以处理含这些元素的体系。如果f电子是局域性很强的(也就是强关联效应),计算出现的问题与一些过渡金属氧化物(比如NiO, V2O3和FeO等)时的一样。
3. POSCAR
描述所计算体系的晶胞参数,原子个数及晶胞中原子的位置,以及分子动力学计算时离子的初始速度(不常用)。
例子:
Si-fcc !注释行,简短描述体系
5.43 !基矢的缩放系数,可认为是晶格常数
0.00 0.50 0.50 !基矢除以缩放系数后的,与上一行的值一起描述基矢
0.50 0.00 0.50
0.50 0.50 0.00
2 !原子个数
Direct !表示原子的坐标是相对于基矢给出的.
0.00 0.00 0.00 !原子的位置
0.25 0.25 0.25
当第七行是C字母开头的,则表示下面的坐标是以卡笛尔坐标系来给出给原子的绝对坐标(被除以了第二行的缩放系数后的坐标值)。比如上面的例子也可以采用下面的方式:
Si-fcc !注释行,简短描述体系
5.43 !基矢的缩放系数,可认为是晶格常数
0.00 0.50 0.50 !基矢除以缩放系数后的,与上一行的值一起描述基矢
0.50 0.00 0.50
0.50 0.50 0.00
2 !原子个数
7Nd 253
Pm 258
Sm 225
Tm 257
Eu 249
Yb 291
Gd 256
Lu 255
U 252 Np 254 Pu 254
Pu_s 211 Ac_s 119 Th_s 169Pa_s 193 U_s 209 Np_s 210
Carti !表示原子的坐标是以卡笛尔坐标系给出的坐标.
0.00 0.00 0.00 !原子的位置
0.25 0.25 0.25
4. KPOINTS
设置布里渊区k点网格取样大小或能带结构计算时沿高对称方向的k点:
a) 手动输入即自定义各个k点的坐标和权重:推荐只在能带计算时用,其他的情况下不采用这种方法。在后面的能带结构计算会详细介绍如何准备手动输入的k点。
例子:
k-points along high symmetry lines !注释行,无特别的意义
11 !沿G-X特殊点之间11个k点
Reciprocal !各k点相对于倒格子基矢来写的
0.00 0.00 0.00 1.00 !k点的坐标和相应的权重
0.05 0.00 0.05 1.00
…….
0.50 0.00 0.50 1.00
b) Line-mode:在计算能带时用(4.6以上版本才支持),不推荐用
例子:
k-points along high symmetry lines !注释行,无特别的意义
10 !沿G-X特殊点之间产生10个k点
Line-mode !程序自动产生特殊k点间的k点
Reciprocal !各k点相对于倒格子基矢来写的
0.00 0.00 0.00 !Gamma
0.50 0.00 0.50 !X
提示: 如果k点是相对于卡笛尔直角坐标系,则第四行改为Cartesian(以字母c开头 的任何词都可以)
c) 程序自动产生k点:最常用的,定义网格取样大小
例子:
Automatic generation !注释行
0 !自动产生k点,这一行必须设为0
Monhkorst-Pack !Monhkorst-Pack方法产生k点
9 9 9 !在各个基矢方向上分割各基矢的点数
0.0 0.0 0.0 !是否移动网格点以及移动多少(这里不移动)
提示:一般各基矢方向上的分割数为奇数,使得产生的k点是以Gamma点为中心的。根据基矢的长短来设置合适的分割数。
针对六角晶系:采用Gamma centered网格
例子:
Automatic generation !注释行
0 !自动产生k点,这一行必须设为0
Gamma !明确定义以Gamma点为中心,根据M-P方法产生k点
9 9 7
8
0.0 0.0 0.0
三、VASP的主要输出文件
VASP的输出文件主要有OUTCAR, CHG, CHGCAR, WAVECAR, DOSCAR, EIGENVAL,
OSZICAR, CONTCAR, PCDAT, IBZKPT, XDATCAR。
1、OUTCAR
OUTCAR文件包含了vasp计算后得到的绝大部分结果,每步迭代的详细情况。下面介绍如何从OUTCAR取出一些有用的信息:
a) 查看所计算体系的体积,使用下面的命令
grep ‘volume’ OUTCAR
得到的结果如下
volume/ion in A,a.u. = 32.92 222.17
volume of cell : 65.84
第一行给出体系的体积分别以Å3/atom, a.u.3/atom为单位给出的。
第二行给出体系的体积是以Å3/unit cell为单位给出的。
b) 查看所计算体系的总能,使用下面的命令
当ISMEAR = -5时,Free energy TOTEN是与energy without entropy是相等,则用
grep ‘TOTEN’ OUTCAR
得到结果如下
free energy TOTEN = -7.910804 eV
当ISMEAR等于其他的值时,Free energy TOTEN是与energy without entropy是不相等,则用
grep ‘entropy=’ OUTCAR
得到结果如下
energy without entropy= -7.910804 energy(sigma->0) = -7.910804
在计算体系的结合能时,体系的总能取为energy without entropy后面的值。
(如何计算体系的结合能,在后面会详细介绍)
c) 查看所计算体系的费米能级,使用下面的命令
grep 'Fermi' OUTCAR | tail -1
得到的结果为
BZINTS: Fermi energy: 6.171330; 20.000000 electrons
9
上一行中第一个数就是体系的费米能级,第二个数就是体系的总价电子数。
注释:对半导体的体系,VASP取价带顶作为费米能级。对呈现金属性的体系,费米能级就是该体系的真实(具有物理意义的)费米能级。
d) 查看所计算体系的倒格子基矢
在采用vi对OUTCAR编辑时,用下面的命令来查找
g/reciprocal lattice vectors 或 g/recip
e) 查看所计算体系中原子的受力情况
在采用vi对OUTCAR编辑时,用下面的命令来查找
g/TOTAL-FORCE
原子所受的力的单位是eV/angstrom。
2、CHG和CHGCAR
这两个都是给出了体系的电荷密度文件,它们的内容是相同的,前者给出的数据的精度要比后者的精度略低一些。下面是CHGCAR文件的例子:
Au-Zn_zig
1.00
15.000000 0.000000 0.000000
0.000000 15.000000 0.000000
0.000000 0.000000 6.600000
1 1
Direct
0.000000 0.000000 0.000000
0.000000 0.000079 0.500000
160 160 72
0.E+05 0.E+05 0.E+05 0.E+05 0.E+05
0.88581841033E+04 0.63620171557E+04 0.42583169365E+04 0.26537018923E+04 0.E+04
...................................
此文件的头9行给出的体系的晶格参数,与POSCAR中的内容基本相同,在11行中三个整数分别是NGX, NGY, NGZ的值,它们表示在三个基矢方向上,对所计算的原胞进行分割,得到NGX * NGY * NGZ个点,所计算原胞中的电荷密度用一个三维矩阵A(NGX, NGY, NGZ)来存贮。
这两个文件在每步迭代过程中都会被更新(除了在INCAR文件中有设置ICHGAR=11或12外)。经过迭代后得到的自洽的CHG和CHGCAR可以用来画图分析面电荷密度分布(如何做,在后面会详细介绍)。在后面步骤中能带结构和态密度时,所读入的电荷密度文件CHG和CHGCAR必须是经过迭代自洽得到的文件。
10
3、DOSCAR和EIGENVAL
DOSCAR给出的所计算体系的电子态密度,EIGENVAL给出的是所计算体系的本征值。这两个文件中的能量值,都是绝对的,而不是以费米能级作为参考零点。
4、其他
WAVECAR给出的是所计算体系的电子波函数,二进制文件,不可编辑。
OSZICAR每次迭代或离子移动情况的简单汇总。
CONTCAR给出的离子进行驰豫时,每次移动后体系的晶格参数,与POSCAR的内容相同。在对体系进行驰豫或分子动力学计算时,最后得到的CONTCAR可以直接拷贝成POSCAR,进行后面的计算。
PCDAT和XDATCAR给出了有关分子动力学模拟中的一些结果,比如pair correlation 函数。
IBZKPT给出的是不可约布里渊区k点的坐标。
lA
参数设置与选择的技巧
下面对一些主要的关键词的设置进行说明。
1、ENCUT
平面波的切断动能。采用默认值还是手动的输入。推荐的做法是采用后者,在任何性质的计算之前,进行ENCUT收敛情况的计算,由此来确定一个合适的切断动能值,然后手动地设置。下面以金刚石结构Si并采用USPP、LDA为例,进行说明:
其POSCAR文件为:
Si-Diamond
5.430
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
0.0 0.0 0.0
0.25 0.25 0.25
其KPOINTS文件:
Automatic generation
0
Monhkorst-Pack
9 9 9
0.0 0.0 0.0
11
用来确定ENCUT的脚本程序为run_ecut,其内容为:
#!/bin/sh
rm WAVECAR
for i in 150 200 250 300 350 400
do
cat > INCAR <
SYSTEM = Si-Diamond
ENCUT = $i
ISTART = 0 ; ICHARG = 2
ISMEAR = -5
PREC = Accurate
!
echo "ENCUT = $i eV" ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f n", $5 }'`
echo $i $E >>comment
done
计算完成后得到comment文件,它列出了在每个ENCUT时计算得到相应的总能。内容如下:
150 -11.900655
200 -11.938864
250 -11.944599
300 -11.945248
350 -11.945503
400 -11.945622
总能变化在0.001eV左右就足够了。在这个例子中,因此,选择ENCUT = 250 eV。
注意循环的第一个值,一般取POTCAR中ENMAX(如果有多个ENMAX,则用其中最大的一个),循环的间隔一般取为50 eV。
另外在对体系的变体积的结构优化时,最好保证ENCUT是ENMAX的1.3倍,以便得到一个合理的精度。但是采用以上方法经过优化得到的ENCUT一般可以满足此条件。
2、PREC
PREC是控制计算精度最重要的一个参数,它决定了ENCUT、FFT网格大小、和ROPT的默认值。可能的取值为Low, Medium, High, Normal, Accurate (后两个只能在4.5以上版本中才起作用)。在一般的计算时推荐:4.5版本中用Normal,4.4版本中用Medium。当要提高力和Stress tensor的计算精度时,可以采用High或Accurate,并手动设置ENCUT的值。
3、 EDIFF和EDIFFG
EDIFF是电子结构部分自洽迭代循环时,判断是否自洽了的条件,上次和当前两次迭代中总能和本征值的变化都小于EDIFF,则电子结构部分迭代循环停止。如果EDIFF = 0,则进行NELM步迭代后停止迭代。默认值为1E-4,一般情况没有必要设置更小的值。
EDIFFG是控制离子部分的驰豫,当离子驰豫在上步和当前步中的总能变化小于EDIFFG,
12
则离子驰豫停止。其默认值为EDIFF*10。注意的是:只有EDIFFG为负数,才是用来控制离子驰豫时,离子或原子所受的力。EDIFFG也可以0,则表示离子驰豫NSW步后就停止。EDIFFG在分子动力学中也不起作用。
4、 ISTART和ICHARG
这两个关键词分别定义了如何构建初始的波函数和电荷密度、读入上一次的波函数和电荷密度。VASP的manual上讲了多种情况,这里推荐的做法是:在进行能带结构、电子态密度等性质的计算时,设置ISTART = 1, ICHARG = 11;其他的情况,一般都设置ISTART = 0,
ICHARG = 2。如果由于断电或其他情况,程序停止运行了,但是又想接着计算,此时在INCAR设置ISTART = 1,ICHARG = 1,其他的参数不变,文件也不用动。
5、GGA和VOSKOWN
GGA关键词表示交换关联函数采用广义梯度近似。当GGA = 91时,表示采用Perdew -Wang
91交换关联函数;当GGA = PE,表示采用Perdew-Burke-Ernzerhof 交换关联函数(只能用在4.5以上版本)。GGA的选择一定要与赝势的类型相一致,也就是说在采用LDA赝势时,不能定义GGA,另外采用PW91的赝势,不能定义GGA = PE等等。VOSKOWN是表示在处理交换关联函数时,采用何种内插公式。当VOSKOWN = 1,采用Vosko Wilk and Nusair提出的内插公式,它一般是来处理GGA中关联函数的,因此在采用PW91的GGA时,应设置VOSKOWN = 1。其他的情况下,可以不必设置VOSKOWN,而由程序采用默认值。
6、ISIF
ISIF是一个非常有用的参数,用来控制结构参数的优化。当IBRION = 0时,其默认值为0,其他情况下为2。ISIF可能的取值以及相应的意思为:
ISIF
0
1
2
3
4
5
6
7
计算离子计算原胞的离子位置所受的力
stress tensor
驰豫
是
是
是
是
是
是
是
是
否
trace only
是
是
是
是
是
是
是
是
是
是
是
否
否
否
改变原胞的形状
否
否
否
是
是
是
是
否
改变原胞的体积
否
否
否
是
否
否
是
是
13
Trace only表示仅有总压力是正确的,总压力也是在OUTCAR文件中这一行“external pressure
= ... kB”给出的。在对原胞的体积或形状进行优化时,ENCUT要略取的大一些(比如取为1.3 * ENCUT的默认值或者设置PREC = High)以消除Pulay Stress导致的误差。
7、ISMEAR和SIGMA
ISMEAR用来确定如何或用何种方法来设置每个波函数的部分占有数fnk。在采用有限温度方法设置fnk时,smearing方法中的smearing宽度σ。它们的默认值分别为ISMEAR = 1,
SIGMA = 0.2。可能的取值为-5, -4, -3, -2, -1, 0, N(N表示正整数),一般很少用-2和-3。
ISMEAR = -5,表示采用Blöchl修正的四面体方法。
ISMEAR = -4,表示采用四面体方法,但是没有Blöchl修正。
ISMEAR = -1,表示采用Fermi-Dirac smearing方法。
ISMEAR = 0,表示采用Gaussian smearing方法。
ISMEAR = N,表示采用Methfessel-Paxton smearing方法,其中N是表示此方法中的阶数。一般情况N取1和2就好,而且大多数情况,N = 1和2给出的结果很接近。
注意:进行任何的静态计算或态密度计算,且k点数目大于4时,取ISMEAR = -5;当由于原胞较大而k点数目较少(小于4个)时,取ISMEAR = 0,并设置一个合适的SIGMA值。另外对半导体或绝缘体的计算(不论是静态还是结构优化),取ISMEAR = -5;当体系呈现金属性时,取ISMEAR = 1和2,以及设置一个合适的SIGMA值。在进行能带结构计算时,ISMEAR和SIGMA用默认值就好。一般说来,无论是对何种体系,进行何种性质的计算,采用ISMEAR = 0,并选择一个合适的SIGMA值都能得到合理的结果。
当采用ISMEAR = 0或N时,如何优化选择SIGMA的值?
以fcc结构Al,采用LDA、USPP为例来进行说明:
其POSCAR文件为
Al-fcc
5.430
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
0.0 0.0 0.0
0.25 0.25 0.25
其KPOINTS文件:
Automatic generation
0
Monhkorst-Pack
9 9 9
14
0.0 0.0 0.0
用来确定SIGMA的脚本程序为run_sigma,其内容为:
#!/bin/sh
rm WAVECAR
for i in 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30
do
cat > INCAR <
SYSTEM = Al-fcc
ENCUT = 250
ISTART = 0 ; ICHARG = 2
ISMEAR = 0; SIGMA = $i
PREC = Accurate
!
echo "SIGMA = $i eV" ; time vasp
TS=`grep "EENTRO" OUTCAR | tail -1 | awk '{printf "%12.6f n", $5 }'`
echo $i $TS >>comment
done
计算完后得到comment文件,其内容为:
0.10 -0.003426
0.12 -0.004408
0.14 -0.005645
0.16 -0.007127
0.18 -0.008833
0.20 -0.010763
0.22 -0.012928
0.24 -0.015336
0.26 -0.017988
0.28 -0.020878
0.30 -0.023999
选择entropy T*S EENTRO值中最小的那个所对应的SIGMA。此例子中,则选择SIGMA =
0.1。当ISMEAR = 1或2时,也可以按这个例子来进行。另外,当k点数目变化后,此SIGMA
值也会要再进行优化。
9、RWIGS
Wigner Seitz半径,用在计算分波态密度以及每根能带对应的波函数按spd和位置投影时。RWIGS = 1.2 1.5 … … 按POSCAR文件中每类原子的顺序相应地给出。尽管在VASP的manual中给出了一个总的原则:调整RWIGS的值,并计算后,检查在OUTCAR文件中,每类原子的Wigner Seitz球的体积之和应略接近与原胞的体积。当体系中有多类原子时,一般很难调整,通常就直接取POTCAR文件中以Å为单位的RWIGS值。
10、k点数目或k-mesh大小的优化
15
以fcc结构Al的计算为例进行说明:
INCAR以一般做静态计算时的情况来设置。
SYSTEM = Al-fcc
ENCUT = 250
ISTART = 0 ; ICHARG = 2
ISMEAR = -5
PREC = Accurate
这个优化的过程可以用下面的脚本程序run_k来完成:
#!/bin/sh
rm WAVECAR
for i in 5 7 9 11 13 15
do
cat > KPOINTS <
Automatic generation
0
Monhkorst-Pack
$i $i $i
0.0 0.0 0.0
!
echo "k mesh = $i x $i x $i" ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f n", $5 }'`
KP=`grep "irreducible" OUTCAR | tail -1 | awk '{printf "%5i n", $2 }'`
echo $i $KP $E >>comment
done
计算完后得到k点数目与能量的对应值,总能变化在0.001eV左右就非常足够了,然后由此来选择合适的k点数目。
五、材料基态性质的计算方法和步骤
在计算前,要明确采用的是何种赝势;平面波切断动能多大;k点网格多小;当采用Gaussian-,Fermi-smearing方法或Methfessel-Paxton smearing方法时,SIGMA多大;计算所选取的精度PREC;采取何种交换关联函数。
另外,在每步计算完后,要学会文件(INCAR, KPOINTS, POSCAR, OUTCAR以及其他的与所计算的性质相关的文件DOSCAR, EIGENVAL)。比如静态计算完后得到自洽的电荷密度,可以建立目录scf,然后把INCAR, KPOINTS, POSCAR, OUTCAR, CHGCAR, CHG保存下来,这可以采用下面的命令来完成:
mkdir scf
tar czvf CHG*
cp INCAR KPOINTS POSCAR OUTCAR scf/.
16
提示:由于CHGCAR的文件比较大,压缩后保存以减少磁盘空间。当要用到时,把解压就可以用了。
比如计算完能带结构,可以建立目录band,然后把INCAR, KPOINTS, POSCAR, OUTCAR,
EIGENVAL, syml文件保存下来,通过下面的命令来完成:
mkdir band
cp INCAR KPOINTS POSCAR OUTCAR EIGENVAL syml band/.
然后进入到目录band下面用pbnd.x程序来处理EIGENVAL。
比如计算电子态密度,可以建立目录dos,然后把INCAR, KPOINTS, POSCAR, OUTCAR,
DOSCAR文件保存下来,通过下面的命令来完成:
mkdir dos
cp INCAR KPOINTS POSCAR OUTCAR DOSCAR dos/.
然后进入目录dos下面用split_dos小程序来处理分割DOSCAR。
1、单个原子的计算
单个原子的计算有两个目的:1) 检验赝势的好坏;2) 对称性被破坏后自旋极化情况下的原子基态能量,对结合能进行修正。
对1)的情况,在VASP的赝势库,由于VASP是商业化的软件,这些元素的赝势都是经过检验过。一般情况下,只要切断动能ENCUT足够大,以及计算单个原子的原胞的晶格常数足够大,得到的能量值应该在1meV~10meV之间,也就是VASP计算得到的单个原子的能量与原子的参考组态时的能量之差。在VASP所计算得到的总能都是扣去了计算原子的参考组态时得到的能量,也就是POTCAR中EATOM的值。
以计算1个Al的情况为例:
KPOINTS的内容为:
Automatic
0
Gamma
1 1 1
0 0 0
POSCAR的内容为:
atom
15.00
1.00000 .00000 .00000
.00000 1.00000 .00000
.00000 .00000 1.00000
1
Direct
0 0 0
INCAR的内容为:
17
SYSTEM = Al: atom
ENCUT = 250.00 eV
NELMDL = 5 !make five delays till charge mixing
ISMEAR = 0; SIGMA=0.1 !use Gaussian smearing method
计算后得到查看OUTCAR文件中的“energy without entropy”之后的能量值。这个值一般要在1meV~10meV之间。
原胞的大小对所有的元素,取15Å是足够的,对某些元素还可以取的更小些。
对2)的情况,还是以计算单个原子Al的为例进行说明:
INCAR文件的内容为:
SYSTEM = Al: atom
ENCUT = 250
ISYM = 0 ! no symmetry
ISPIN = 2 ! allow for spin polarisation
VOSKOWN = 1 ! this is important, in particular for GGA
ISMEAR = 0 ! Gaussian smearing, otherwise negative occupancies
SIGMA = 0.1 ! intermid. smearing width
AMIX = 0.2 ! mixing set manually
BMIX = 0.0001
NELM = 20 ! 20 electronic steps
ICHARG = 1
连续计算两次,查看OUTCAR文件中的“energy without entropy”之后的能量值,这个值就是用来修正体材料的结合能的。
原胞的大小,与1)情况中的相同。上面INCAR中的内容从第3行起后面的设置,可以用在计算其他原子的情况中。
2、结构参数(晶格常数和原子位置参数)的优化
根据要优化的晶胞参数的复杂性可以分为以下几类:
1) 简单的情况:只要优化一个参数即晶格常数a,其步骤如下(以计算fcc结构Al的晶格常数为例进行说明):
a) 准备好INCAR,即定义ENCUT,ISTART = 0,ICHARG = 2,ISMEAR = -5
SYSTEM = Al-fcc
ENCUT = 250
ISTART = 0; ICHARG = 2
ISMEAR = -5
PREC = Accurate
b) 准备好KPOINTS,POTCAR(为USPP, LDA)
Automatic generation
0
Monhkorst-Pack
9 9 9
0.0 0.0 0.0
18
c) 准备好POSCAR文件,以晶格常数实验值aexp为基础,在aexp左右计算10个点得到Volume-Etotal的数据。这个可以通过脚本程序run_a0来完成
#!/bin/sh
rm WAVECAR
for i in 3.80 3.85 3.90 3.95 4.00 4.05 4.10 4.15 4.20 4.25 4.30
do
cat > POSCAR <
Al-fcc
$i
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1
Direct
0.0 0.0 0.0
!
echo “ a = $i angstrom “; time vasp
E=`grep “TOTEN” OUTCAR | tail -1 | awk ‘{printf “%12.6f n”, $5 }’`
V=`grep “volume” OUTCAR | tail -1 | awk ‘{printf “%12.4f n” , $5}’`
echo $V $E >>
done
得到的文件,其内容如下:
13.7200 -4.094976
14.2700 -4.137590
14.8300 -4.163643
15.4100 -4.176403
16.0000 -4.176731
16.6100 -4.166067
17.2300 -4.145942
17.8700 -4.117937
18.5200 -4.083448
19.1900 -4.043300
19.8800 -3.998039
其中第一列数据是体积,单位为Å3,第二列数据是能量,单位为eV。
d) 采用Rose公式或Birch-Murnaghan状态方程拟合得到晶格常数。
2) 复杂的情况:含两个以上的参数,比如四角或六角晶系(a,c),正交晶系(a,b,c);以及含有原子位置参数需要优化,步骤为(以计算六角结构Mg的晶格参数为例进行说明):
a) 以实验的晶格结构参数为基础,做好POSCAR,先确定好ENCUT,k-mesh大小,SIGMA的值,准备一个名为的文件,其内容大致如下:
SYSTEM = Mg-hex
ENCUT = 250
ISTART = 0; ICHARG = 2
19
ISMEAR = 1; SIGMA = 0.2
NSW = 60; IBRION = 2
ISIF = 5
POTIM = 0.2
EDIFF = 1E-5; EDIFFG = -1E-3
PREC = Accurate
再准备一个名为的文件,其内容大致如下:
SYSTEM = Mg-hex
ENCUT = 250
ISTART = 0; ICHARG = 2
ISMEAR = -5
PREC = Accurate
其中POSCAR的内容如下:
Auto generation
0
Gamma
9 9 7
0.0 0.0 0.0
b) 先进行一次体积保持不变的离子驰豫的计算(通过ISIF来设置,此时ISIF可能的取值为2,4或5)
ISIF的选择根据所要优化的结构参数的来进行选择,见上一部分对ISIF的说明。其中“改变原胞的形状”,也就是调整原胞中c/a和b/a的值。
c) 再把优化得到的CONTCAR拷贝成POSCAR,进行一次静态的计算
d) 对a的值取10个左右的点,每个点重复上面两步,得到静态计算下的Volume-Etot关系。这三步可以通过运行脚本程序run_cell来进行,其中run_cell的内容如下:
#!/bin/sh
rm WAVECAR
for i in 2.81 2.91 3.01 3.11 3.21 3.31 3.41 3.51 3.61 3.71
do
cat > POSCAR <
Mg-hex
$i
0.0 -1.0 0.0
0.866 0.5 0.0
0.0 0.0 1.6230529595
2
Direct
0.6666666666666667 0.3333333333333333 0.750
0.3333333333333333 0.6666666666666667 0.250
!
cp INCAR
echo "a = $i angstrom " ; time vasp
cp CONTCAR POSCAR
20
cp INCAR
echo "a = $i angstrom " ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f n", $5 }'`
V=`grep "volume" OUTCAR | tail -1 | awk '{printf "%12.4f n", $5 }'`
echo $V $E >>
done
在run_cell运行完后,得到文件。
e) 采用状态方程拟合得到平衡状态下的体积,体弹性模量
f) 在该体积下,重复上面b)和c)步,得到平衡状态下的其他晶胞参数。这一步也就是:在得到了E(V)曲线后,通过状态方程拟合得到平衡状态下的体积,计算出上面脚本中变量$i的值,并改变$i的循环值,再运行run_cell计算一次,得到其他的结构参数c和位置u.。
另外一种对体系的结构参数进行一次性型的计算(这种方法一般是用来估计的,计算得到较合理,但是精度不高)。这通过设置ISIF来进行的。还是以计算六角结构Mg为例:
计算时的INCAR文件为:
SYSTEM = Mg-hex
ENCUT = 250
ISTART = 0; ICHARG = 2
ISMEAR = 1; SIGMA = 0.2
NSW = 60; IBRION = 2
ISIF = 3
POTIM = 0.2
EDIFF = 1E-6; EDIFFG = -1E-3
PREC = Accurate
注释:此时可以把EDIFF和EDIFFG的精度提高一些以得到更准确的晶格参数。
KPOINTS与前面的相同。POSCAR的内容为:
Mg-hex
3.21
0.0 -1.0 0.0
0.866 0.5 0.0
0.0 0.0 1.6230529595
2
Direct
0.6666666666666667 0.3333333333333333 0.750
0.3333333333333333 0.6666666666666667 0.250
最后计算完后,得到的CONTCAR文件就包含优化后的晶格参数。
这样也可以比较采用这两种方法得到的晶格参数究竟差多少。
3、结合能
VASP计算得到的总能已经减去了在以原子参考组态计算得到的原子能量(也就是构造赝势
21
时,得到的总能,对应于POTCAR文件中的EATOM)。要得到准确的结合能,还需减去前面单个原子计算中的第2)种情况计算得到的修正值。
4、自洽的电荷密度
再优化得到了晶胞参数后,进行静态的计算就可以得到自洽的电荷密度,并要保存下来,在后面计算其他的性质时要用到;另外也可以根据它画出面电荷密度图,分析原子间的键合作用。步骤为(并以计算fcc结构Al为例进行说明):
a) 准备好INCAR,即定义ENCUT,ISTART=0,ICHARG=2,ISMEAR=-5
SYSTEM = Al-fcc
ENCUT = 250
ISTART = 0; ICHARG = 2
ISMEAR = -5
PREC = Accurate
b) 准备好KPOINTS和POTCAR
Automatic generation
0
Monhkorst-Pack
9 9 9
0.0 0.0 0.0
(这个是KPOINTS文件中的内容)
c) 准备好POSCAR文件或以优化的晶格参数作为基础,把优化得到的CONTCAR拷贝成POSCAR。
Al-fcc
3.975
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1
Direct
0.0 0.0 0.0
d) 提交运行:用命令nohup time vasp&
e) 当计算完成后,保存CHGCAR和CHG:用命令tar czvf CHG*
f) 用命令cp CHGCAR ,并仅在文件中第一行后加入_P1_charge / x x
x … /,/ x x x …/按POSCAR文件中每类原子的名称给写出。然后用VENUS软件打开文件,进行面电荷密度的分析。
Al-fcc_P1_charge / Al /
3.975
0.000000 0.500000 0.500000
0.500000 0.000000 0.500000
0.500000 0.500000 0.000000
22
1
Direct
0.000000 0.000000 0.000000
28 28 28
………………
5、能带结构
计算材料的能带结构即色散曲线E(k),步骤为(并以计算fcc结构Al的能带结构为例进行说明):
a) 根据特殊k点的走向,选取特殊k点及特殊k点间的分割点数,准备好产生k点的输入文件syml
6 !特殊k点的个数
20 20 20 10 20 !特殊k点间的分割点数
X 0.5 0.0 0.5 !特殊k点的坐标,相对于倒格子矢量
G 0.0 0.0 0.0
L 0.5 0.5 0.5
W 0.5 0.25 0.75
K 0.375 0.375 0.75
G 0.0 0.0 0.0 !下面三行,前三列是正格子基矢,后三列是倒格子基矢
0.000000000 1.987500000 1.987500000 -0.251572327 0.251572327 0.251572327
1.987500000 0.000000000 1.987500000 0.251572327 -0.251572327 0.251572327
1.987500000 1.987500000 0.000000000 0.251572327 0.251572327 -0.251572327
-20.0 15.0 !在画能带结构时,每个特殊k点所对应的竖线的能量范围
7.068339 !费米能级
b) 用程序gk.x产生k点,得到KPOINTS文件
注释:程序gk.x是由gk.f文件编译后得到的目标文件,其输入文件为syml,输出文件为KPOINTS, 。
c)紧接着利用前面计算得到的自洽电荷密度作一次非自洽的计算
采用命令解压保存的电荷密度文件:tar xzvf
另外设置ISTART=1, ICHARG=11, 并增加NBANDS的值,ISMEAR采用默认值
SYSTEM = Al-fcc
ENCUT = 250
ISTART = 1; ICHARG = 11
#ISMEAR = -5
NBANDS = 12
PREC = Accurate
计算完后得到本征值文件EIGENVAL。
注意:对于4.4系列版本,在计算能带结构时设置NBANDS的值应该与计算自洽的电荷密度时设置的NBADS一致。对4.5以上版本,可以不一致。
d) 从自洽电荷密度计算得到的OUTCAR文件中找到倒格子矢量和费米能级,并粘贴到syml
23
文件中,然后用程序pbnd.x把EIGENVAL转换为成(本征值,并以费米能级为参考零点)和(用来画竖线),然后用软件origin画图。
注释:程序pbnf.x是通过编译pbnd.f得到的可执行文件,其输入文件为EIGENVAL和syml,输出文件为BANDS、和。pbnd.f可以处理自旋极化情况下计算得到的EIGENVAL,不再输出而是和这两个文件,分别对自旋向上和向下的能带。
提示:在计算能带结构时,采用ISMEAR = 0或1对结果的影响非常小,可以认为是一样的。但是不能采用ISMEAR = -5 或-4。
6、电子态密度
计算材料的电子态密度可以包括总态密度和分波态密度,步骤为(以计算fcc结构Al的态密度为例子进行说明):
a) 准备好KPOINTS文件,增加k点网格
Automatic generation
0
Mohkorst-Pack
19 19 19
0.0 0.0 0.0
b) 从POTCAR文件中找到各类原子的RWIGS(vi编辑POTCAR,并用命令g/ RWIGS)
c) 准备好INCAR文件(设置ISTART=1, ICHARG=11, ISMEAR=-5以及RWIGS)
SYSTEM = Al-fcc
ENCUT = 250
ISTART = 1; ICHARG = 11
ISMEAR = -5
RWIGS = 1.402
PREC = Accurate
d) 利用前面计算得到的自洽电荷密度,进行一次非自洽计算
tar xzvf
nohup time vasp&
计算完后,得到包含了态密度值的DOSCAR文件,
e) 采用split_dos对态密度文件DOSCAR进行分割,得到总态密度DOS0,各个原子的分波态密度DOS1,DOS2……。另外在运行split_dos程序对DOSCAR文件分割时,要保证当前目录下有对应的OUTCAR和POSCAR文件。
分割后的DOS0,DOS1…等文件的能量值是以费米能级作为能量参考零点。DOS0的第一列数据是能量值,单位为eV;第二列数据是总态密度的值,单位State/ cell;第三列数据是总态密度的积分值,也就是电子数,单位为electrons。DOS1是第一个原子的分波态密度值,其中的第一列数据是能量值,单位为eV;第二、三、四列数据分别对应于s、p、d
24
态的分波态密度值,单位为State/。其他的DOS文件与DOS1类似。
六、材料磁性性质的计算
磁性的计算,其实与非磁性的计算相比,就只需在INCAR中加入ISPIN = 2以及设置各类原子的初始磁矩,这通过MAGMOM来设置。更复杂的磁性性质的计算,包括noncollinear磁性、spin orbit相互作用和Spin sprial磁性,需要再增加其他的关键词。下面主要讲的是如何进行简单的磁性计算。根据设置MAGMOM的不同来确定计算材料的铁磁、反铁磁以及亚铁磁性质。以计算fcc结构Ni的铁磁性为例进行说明(在例子,采用的是PBE、PAW势):
提示:作磁性计算时,强烈推荐采用PAW势,得到的结果会更准确些。
其晶格参数、基态性质的计算基本与非磁性时的计算相同,只需在INCAR文件中加入SPIN
= 2以及设置MAGMOM的值。
INCAR文件的内容为:
SYSTEM = Ni-FM
ISTART = 0; ICHARG = 2
ENCUT = 350 eV
ISMEAR = -5
GGA = PE; VOSKOWN = 1
ISPIN = 2
MAGMOM = 1
PREC = Accurate
注释:当采用GGA = 91时,强烈推荐VOSKOWN = 1加上。当采用GGA = PE,VOSKOWN
= 1,可加可不加,此时如果没有加上VOSKOWN = 1,程序默认也是采用Vosko Wilk and
Nusair提出的内插公式来处理关联部分。
MAGMOM的设置要对应于POSCAR文件中每类原子。进行铁磁性质的计算时,MAGMOM要设置成相同的值,在INCAR文件中,也可以不设置,程序会默认设置为1。
KPOINTS和POSCAR与进行非磁性时的一样。
Automatic generation
0
Monhkorst-Pack
9 9 9
0.0 0.0 0.0
(这个是KPOINTS文件中的内容)
Ni-FM
3.52
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
25
1
Direct
0.0 0.0 0.0
(这个是POSCAR文件的内容)
计算完,从OSZICAR最后一行可以找到体系总的磁矩是0.567µB。这个也可以从通过grep
‘magnetization‘ OUTCAR在OUTCAR文件中找到(看“number of electron …
magnetization …”这一行的数据)。
在后面计算能带结构和电子态密度时,分别会得到自旋向上和自旋向下的两部分。
七、表面体系的计算
在作表面的性质计算时,现在一般都是采用slab(“晶层”或“薄片”)模型来模拟表面体系的。因此表面体系的计算大致可以分为四个大的步骤:1、材料体性质的计算;2、slab模型的构造;3、表面体系的结构优化;4、表面体系性质的计算。
1、材料体性质的计算
这一步包含了前面材料基态性质的计算。主要是为了确定后面在进行表面计算时所需要的一些参数:ENCUT,当采用ISMEAR = 1或0时的 SIGMA,以及体材料的晶格参数(因为构造slab模型是以体材料的晶格参数作为基础来进行的)。如何确定这些参数,可以参考前面所介绍的方法和步骤。其中SIGMA的优化是必须的,因为后面对表面体系的结构进行优化时,smearing方法一般都是采用Gaussian方法或Methfessel-Paxton smearing方法。
2、slab模型的构造
在知道了体材料的晶格参数,以及明确了要模拟什么样的表面(也就是表面的密勒指数以及表面的二维周期性),就可以开始构造该表面的slab模型了。在构造slab模型时,还有两个重要的参数,就是真空层(Vacuum layer)和原子层的厚度,这是因为slab模型就是由原子层和真空层所组成的。真空层和原子层要取多厚,这要经过试用不同的厚度得到,看它们对总能的影响,然后选择合适的厚度。
3、表面体系的结构优化
在对表面体系的结构进行优化前,还需要对k点数目或k mesh大小进行优化,这个的优化也可以参考前面的介绍。在对表面体系的结构进行优化时,主要是对原子的位置进行优化,而不再对超原胞(Slab模型得到的)大小进行优化,一般采用的是Selective Dynamic(也就是有选择性的位置驰豫)。这是在POSCAR进行设置的,另外在INCAR文件也应加入控制离子驰豫的关键词。下面以优化Al(100)-p(1x1)为例进行说明:
INCAR文件的内容为:
26
SYSTEM = Al(100)-p(1x1)
ENCUT = 200
ISMEAR = 1; SIGMA = 0.20
ISTART = 0; ICHARG = 2
EDIFF = 1E-5; EDIFFG = -1.0E-3
NSW = 60; IBRION = 2
POTIM = 0.1
PREC= Accurate
KPOINTS文件的内容为:
auto
0
Monkhorst-Pack
1 11 11
0.0 0.0 0.0
POSCAR文件的内容为:
Al(100)-p(1x1)
1.00000
0. 2. -2.
0. 2. 2.
22.1485000000 0. 0.
7
Selective dynamics
Direct
0. 0. 0. F F T
0. 0. 0.1828340520 F F F
0. 0. 0.3656681039 F F F
0. 0. 0.5485021559 F F T
0.5000000000 0.5000000000 0. F F T
0.5000000000 0.5000000000 0.2742510780 F F F
0.5000000000 0.5000000000 0.4570851299 F F T
其中对原子层上下各两层原子进行驰豫,中间三层原子位置固定。
在INCAR文件中的EDIFF,EDIFFG,NSW控制原子位置驰豫的步数。当原子所受的力小于EDIFFG时,原子位置就停止移动。得到的CONTCAR文件,就是驰豫得到的最后位置。
原子受力的情况,可以在OUTCAR文件中查找TOTAL-FORCE来查看。
4、表面体系性质的计算。
在得到了优化的结构,就可以进行一系列性质的计算,其步骤与前面体材料性质的计算一样。
提示:无论是对体材料还是表面体系的结构优化,在结构优化完后,还需进行进行静态的计算,以得到自洽的电荷密度,再进行后面的性质计算。结构优化完得到的电荷密度文件不可用在后面的性质计算中。在计算功函数和进行STM模拟时,需另加其他的关键词。如对此有兴趣,后面将作专题介绍。
27
”ªA tools中小程序的说明
此部分是对tools中的一些小程序进行说明:
1、murn.f
这个程序是采用Murn状态方程来拟合晶格常数和计算体弹性模量的。从internet网上收集到的。
编译:
g77 -o murn.x murn.f
得到可执行文件murn.x。
使用:
其输入文件为inp.m,inp.m的内容以及格式为:
2
0.25
6.00 7.45 50
8
6.0849 -.62120891E+01
6.2739 -.64553428E+01
6.4629 -.65246916E+01
6.6518 -.64694519E+01
6.8408 -.63284020E+01
7.0298 -.61271095E+01
7.2188 -.58843994E+01
7.4077 -.56183090E+01
-----------------------------------------
第一行表示下面输入的能量的单位是采用eV、Ry或Hartree。当采用的是Ry,则为1;如果是eV,则用2;如果是Hartree,则用3。
第二行可以看成是原胞的体积与晶胞的体积之比。这里晶胞的体积计算公式为:晶格常数的三次方,原胞的体积计算公式按固体物理教科书中的方法来计算。比如所计算的体系是fcc,则为0.25;如果是bcc,则为0.5;如果是六角的sqrt(3.0)/2 * c/a。其他的自己去推算了。
第三行是用来拟合的晶格常数时晶格常数的范围,以及点数。第一个数是晶格常数的最小值,第二个数是晶格常数的最大值,第三个是拟合多少个点,这个数一般取30~50。
第四行是所计算了多少个晶格常数的能量值用来拟合。
28
从第五行起,是计算得到的晶格常数与能量的一一对应值,第一列是晶格常数;第二列是能量值,它的单位要与前面第一行的确定的单位一致。晶格常数-能量值的对数也要与第四行的数一致。
文件准备好之后,使用的方法为:murn.x
使用murn.x拟合后得到的拟合值写在文件中,拟合的情况及重要的结果写到out.m文件中。
用vi编辑out.m,然后用命令g/alat=来查找拟合得到的晶格常数(单位为a.u.)和体弹性模量(单位为Mbar)。
2、gk.f和pbnd.f
gk.f是用来产生计算能带结构时所需要的k点,其输入文件为syml,在手册中有详细介绍每行的意思,这里就不再赘述。输出文件为KPOINTS和。说明一下,在syml中特殊k点的总数不能超过10个。如遇到超出10,则可以把gk.f中有关定义特殊k点的数组的维数调大。
编译:
g77 -o gk.x gk.f
pbnd.f是用来把本征值文件EIGENVAL转换为可以用origin软件来画图的数据。其输入文件为syml和EIGENVAL。输出文件为(或和)和。说明一下,pbnd.f能写出的能带数,默认最大是100个,如果超过了100,可以在pbnd.f中调大定义的值。
编译:
g77 -o pbnd.x pbnd.f
另外还附带针对fcc、bcc、sc和六角晶系的syml,分别名为, , ,
。使用时,把相应的拷贝成syml。
3、split_dos和vp
这两个要一起用的,都是csh脚本程序,其中运行split_dos要调用vp。从internet网上收集的。
它是用来分割DOSCAR,把DOSCAR分解成每个原子的,以方便用origin来画图。其使用的说明也可以参见前面的介绍。
使用时,在你的当前主目录建立一个bin的目录。比如你的用户名为xxxx,则在/home/xxxx目录下,建立bin目录(mkdir /home/xxxx/bin),然后把split_dos和vp放到该目录下。
29


发布评论