VASP的个人经验手册

使用VASP的个人经验手册

侯柱锋

厦门大学物理系2004届博士 E-mail: zhufhou@yahoo.com

2004/06/22

本手册纯属个人使用VASP后的心得和经验总结，版权属于本手册的作者及厦门大学物理系计算物理实验室（Group leader: 朱梓忠教授）。未经许可，不准在网上传阅。文中提到的一些小程序，可以提供使用。在参考的过程，如遇到不清楚或含糊的地方，可以参考VASP的英文manual或email给我。如认为本手册某些地方需要更正或修改的，请email给我。当在使用VASP的过程中遇到问题，也可以email给我，大家一起学习VASP的使用，挖掘和掌握VASP强大的功能。本手册参考了VASP的英文manual、G.Kresse的报告以及从internet网上收集的资料。

本手册大致有以下几个内容： A

程序的编译

“?AVASP的主要输入文件 OAVASP的主要输出文件 lA参数设置与选择的技巧 A

材料基态性质的计算方法和步骤

ZA材料磁性性质的计算 μA表面体系的计算 ”aAtools中小程序的说明 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源代码：vasp.4.4.5.tar.gz和vasp.4.lib.tar.gz 数学库：LAPACK和BLAS (http://www.netlib.org/)， 或mkl（配合intel的fotran编译器用）， 或ATLAS (http://math-atlas.sourceforge.net/)

或Lib GOTO (http://www.cs.utexas.edu/users/flame/goto/)

Fortran编译器：PGI fortran 至少4.0以上版本(http://www.pgroup.com/)，

或Intel的 ifc (8.0以上版本是ifort, http://www.intel.com/software/products/compilers/flin/)，前者可以从网站上下载到15天的试用版本，后者可以从网站下载到免费的版本。或者在国内的个人ftp服务器上搜索它们的破解版本。 本实验室的都有这些软件的备份。

2、下面采用PGI fortan编译器pgf90、ATLAS数学库对VASP4.4.5进行编译

这里假定已经安装好了fortran编译器，所有文件都放在/home/houzf/VASP_SRC目录下，机器的操作系统是Linux: Redhat9.0。

a) 从http://math-atlas.sourceforge.net/下载atlas3.6.0_Linux_P4SSE2.tar.gz，并用如下命令解压：tar xzvf atlas3.6.0_Linux_P4SSE2.tar.gz

解压后得到一个目录Linux_P4SSE2，在此目录下有个lib子目录，该lib子目录中的文件为libatlas.a, libcblas.a, libf77blas.a, liblapack.a, 这些就是编译vasp时所需要的数学库文件之一。 b) 用如下命令解压vasp.4.4.5.tar.gz和vasp.4.lib.tar.gz： tar xzvf vasp.4.4.5.tar.gz tar xzvf vasp.4.lib.tar.gz

解压后分别得到目录vasp.4.4和vasp.4.lib，目录vasp.4.4中文件是vasp的主要源代码，vasp.4.lib是编译vasp时需要的一些特定的数学库程序，在这两个目录中都有编译时所用的makefile文件，针对机器和fortran编译器，选择相应的makefile。

c) 进入vasp.4.lib目录，选择makefile.linux_pg，并把它拷贝成makefile，然后键入make命令开始编译。整个命令如下： cd vasp.4.lib

cp makefile.linux_pg makefile make

编译成功后，得到libdmy.a文件。

d) 退出vasp.4.lib目录，进入vasp.4.4目录，选择makefile.linux_pg，并把它拷贝成makefile，编辑makefile文件，通过修改LIB变量的赋值而采用基于ATLAS的数学库文件，修改的地

2

方和方法是：

在第87和88行前加上#，把这两行注释掉，然后去掉第91,92和93行前的#。 修改前和后的内容为分别为：

LIB = -L../vasp.4.lib -ldmy ../vasp.4.lib/linpack_double.o \\ ../vasp.4.lib/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../vasp.4.lib -ldmy ../vasp.4.lib/linpack_double.o \\ # ../vasp.4.lib/lapack_double.o -L/usr/local/lib \\

# -L$(HOME)/archives/BLAS_OPT/ATLAS/lib/Linux_ATHLONTB/ -lf77blas –latlas #LIB = -L../vasp.4.lib -ldmy ../vasp.4.lib/linpack_double.o \\ # ../vasp.4.lib/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../vasp.4.lib -ldmy ../vasp.4.lib/linpack_double.o \\ ../vasp.4.lib/lapack_double.o -L/usr/local/lib \\ -L../Linux_P4SSE2/lib/ -lf77blas -latlas 修改后保存makefile文件，键入make命令开始编译vasp。整个命令为： cd .. cd vasp.4.4

cp makefile.linux_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 250

O 400

F_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 98

Sb 172

As 208

Te 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 !原子个数

7

Nd 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.00000000000000

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.18441120499E+05 0.17909524567E+05 0.16406959292E+05 0.14179806898E+05 0.11554638997E+05 0.88581841033E+04 0.63620171557E+04 0.42583169365E+04 0.26537018923E+04 0.15676950926E+04

...................................

此文件的头9行给出的体系的晶格参数，与POSCAR中的内容基本相同，在11行中三个整数分别是NGX, NGY, NGZ的值，它们表示在三个基矢方向上，对所计算的原胞进行分割，得到NGX * NGY * NGZ个点，所计算原胞中的电荷密度用一个三维矩阵A(NGX, NGY, NGZ)来存贮。

这两个文件在每步迭代过程中都会被更新（除了在INCAR文件中有设置ICHGAR=11或12外）。经过迭代后得到的自洽的CHG和CHGCAR可以用来画图分析面电荷密度分布（如何做，在后面会详细介绍）。在后面步骤中能带结构和态密度时，所读入的电荷密度文件CHG和CHGCAR必须是经过迭代自洽得到的文件。

10

以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 <

Monhkorst-Pack $i $i $i 0.0 0.0 0.0 !

echo \ time vasp

E=`grep \ '{printf \ $5 }'` KP=`grep \ '{printf \ $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.tgz CHG*

cp INCAR KPOINTS POSCAR OUTCAR chg.tgz scf/.

16

提示：由于CHGCAR的文件比较大，压缩后保存以减少磁盘空间。当要用到时，把chg.tgz解压就可以用了。

比如计算完能带结构，可以建立目录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 <

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 “.6f \\n”, $5 }’` V=`grep “volume” OUTCAR | tail -1 | awk ‘{printf “.4f \\n” , $5}’` echo $V $E >>EtVo.dat done

得到的EtVo.dat文件，其内容如下： 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的值，准备一个名为INCAR.relax的文件，其内容大致如下： 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

再准备一个名为INCAR.static的文件，其内容大致如下： 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 <

0.0 -1.0 0.0 0.8660254037844 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.relax INCAR

echo \ time vasp cp CONTCAR POSCAR

20

cp INCAR.static INCAR

echo \ time vasp

E=`grep \ '{printf \ $5 }'` V=`grep \ \echo $V $E >>EtVo.dat done

在run_cell运行完后，得到EtVo.dat文件。

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.8660254037844 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.tgz CHG*

f) 用命令cp CHGCAR rho.vasp，并仅在rho.vasp文件中第一行后加入_P1_charge / x x x … /，/ x x x …/按POSCAR文件中每类原子的名称给写出。然后用VENUS软件打开rho.vasp文件，进行面电荷密度的分析。 Al-fcc_P1_charge / Al / 3.97500000000000

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, inp.kpt。

c）紧接着利用前面计算得到的自洽电荷密度作一次非自洽的计算 采用命令解压保存的电荷密度文件chg.tgz：tar xzvf chg.tgz

另外设置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转换为成bnd.dat（本征值，并以费米能级为参考零点）和highk.dat（用来画竖线），然后用软件origin画图。

注释：程序pbnf.x是通过编译pbnd.f得到的可执行文件，其输入文件为EIGENVAL和syml，输出文件为BANDS、bnd.dat和highk.dat。pbnd.f可以处理自旋极化情况下计算得到的EIGENVAL，不再输出bnd.dat而是upbnd.dat和dnbnd.dat这两个文件，分别对自旋向上和向下的能带。

提示：在计算能带结构时，采用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 chg.tbz nohup time vasp&

计算完后，得到包含了态密度值的DOSCAR文件，

e) 采用split_dos对态密度文件DOSCAR进行分割，得到总态密度DOS0，各个原子的分波态密度DOS1，DOS2……。另外在运行split_dos程序对DOSCAR文件分割时，要保证当前目录下有对应的OUTCAR和POSCAR文件。

分割后的DOS0，DOS1…等文件的能量值是以费米能级作为能量参考零点。DOS0的第一列数据是能量值，单位为eV；第二列数据是总态密度的值，单位State/eV.unit cell；第三列数据是总态密度的积分值，也就是电子数，单位为electrons。DOS1是第一个原子的分波态密度值，其中的第一列数据是能量值，单位为eV；第二、三、四列数据分别对应于s、p、d

24

态的分波态密度值，单位为State/eV.atom。其他的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.0000000000 2.0247500000 -2.0247500000 0.0000000000 2.0247500000 2.0247500000 22.1485000000 0.0000000000 0.0000000000 7

Selective dynamics Direct

0.0000000000 0.0000000000 0.0000000000 F F T 0.0000000000 0.0000000000 0.1828340520 F F F 0.0000000000 0.0000000000 0.3656681039 F F F 0.0000000000 0.0000000000 0.5485021559 F F T 0.5000000000 0.5000000000 0.0914170260 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

”aA 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 out.m

使用murn.x拟合后得到的拟合值写在文件fout.dat中，拟合的情况及重要的结果写到out.m文件中。

用vi编辑out.m，然后用命令g/alat=来查找拟合得到的晶格常数（单位为a.u.）和体弹性模量（单位为Mbar）。

2、gk.f和pbnd.f

gk.f是用来产生计算能带结构时所需要的k点，其输入文件为syml，在手册中有详细介绍每行的意思，这里就不再赘述。输出文件为KPOINTS和inp.kpt。说明一下，在syml中特殊k点的总数不能超过10个。如遇到超出10，则可以把gk.f中有关定义特殊k点的数组的维数调大。 编译：

g77 -o gk.x gk.f

pbnd.f是用来把本征值文件EIGENVAL转换为可以用origin软件来画图的数据。其输入文件为syml和EIGENVAL。输出文件为bnds.dat(或upbnd.dat和dnbnd.dat)和highk.dat。说明一下，pbnd.f能写出的能带数，默认最大是100个，如果超过了100，可以在pbnd.f中调大定义的值。 编译：

g77 -o pbnd.x pbnd.f

另外还附带针对fcc、bcc、sc和六角晶系的syml，分别名为syml.fcc, syml.bcc, syml.sc, syml.hex。使用时，把相应的拷贝成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

本文来源：https://www.bwwdw.com/article/0777.html