VASP几个计算实例.doc

用VASP计算原子的能量 氢原子的能量为 。在这一节中，我们用VASP计算H原子的能量。对于原子计算，我们可以采用如下的INCAR文件 PREC=ACCURATE NELMDL = 5 make five delays till charge mixing ISMEAR = 0; SIGMA=0.05 use smearing method采用如下的KPOINTS文件。由于增加K点的数目只能改进描述原子间的相互作用，而在单原子计算中并不需要。所以我们只需要一个K点。 Monkhorst Pack 0 Monkhorst Pack 1 1 1 0 0 0采用如下的POSCAR文件 atom 1 15.00000 .00000 .00000 .00000 15.00000 .00000 .00000 .00000 15.00000 1cart 0 0 0采用标准的H的POTCAR 得到结果如下： k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation 1 -6.3145 1.00000 2 -0.0527 0.00000 3 0.4829 0.00000 4 0.4829 0.00000我们可以看到，电子的能级不为 。 Free energy of the ion-electron system (eV) - alpha Z PSCENC = 0.00060791 Ewald energy TEWEN = -1.36188267 -1/2 Hartree DENC = -6.27429270 -V(xc)+E(xc) XCENC = 1.90099128 PAW double counting = 0.00000000 0.00000000 entropy T*S EENTRO = -0.02820948 eigenvalues EBANDS = -6.31447362 atomic energy EATOM = 12.04670449 - free energy TOTEN = -0.03055478 eV energy without entropy = -0.00234530 energy(sigma-0) = -0.01645004我们可以看到 也不等于 。 在上面的计算中有个问题，就是H原子有spin,而在上面的计算中我们并没有考虑到spin。所以如果我们改用LSDA近似，在INCAR中用ISPIN=2的tag，则得到如下结果： k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation 1 -7.2736 1.00000 2 -0.1229 0.00000 3 0.4562 0.00000 4 0.4562 0.00000 5 0.4562 0.00000 spin component 2 k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation 1 -2.4140 0.00000 2 -0.0701 0.00000 3 0.5179 0.00000 4 0.5179 0.00000 5 0.5179 0.00000 Free energy of the ion-electron system (eV) - alpha Z PSCENC = 0.00060791 Ewald energy TEWEN = -1.36188267 -1/2 Hartree DENC = -6.68322940 -V(xc)+E(xc) XCENC = 2.38615430 PAW double counting = 0.00000000 0.00000000 entropy T*S EENTRO = 0.00000000 eigenvalues EBANDS = -7.27361676 atomic energy EATOM = 12.04670449 - free energy TOTEN = -0.88526212 eV energy without entropy = -0.88526212 energy(sigma-0) = -0.88526212氢原子的能量约等于 。可以看到在LDA中如果限制自旋，使能级大概提高了 。但是如何理解所得到的能级，由于用到了赝势，本人并不很清楚如何解释能级意义。 用VASP计算Pd金属的晶格常数 Pd金属的实验上的晶格常数为 。在这里，我们用VASP计算它的晶格常数。首先将Pd所对应的POTCAR文件拷贝到目录下。然后准备好INCAR和KPOINTS文件。POSCAR文件我们将通过一个tcsh的script来产生。 KPOINTS文件可以如下： Monkhorst Pack 0 Monkhorst Pack 11 11 11 0 0 0INCAR文件可以如下： SYSTEM = Pd bulk calculation Startparameter for this run: PREC = Accurate ISTART = 0 job : 0-new 1-cont 2-samecut ICHARG = 2 charge: 1-file 2-atom 10-const ISPIN = 1 spin polarized calculation? Electronic Relaxation 1 EDIFF = 0.1E-03 stopping-criterion for ELM LREAL = .FALSE. real-space projection Ionic relaxation EDIFFG = 0.1E-02 stopping-criterion for IOM NSW = 0 number of steps for IOM IBRION = 2 ionic relax: 0-MD 1-quasi-New 2-CG ISIF = 2 stress and relaxation POTIM = 0.10 time-step for ionic-motion TEIN = 0.0 initial temperature TEBEG = 0.0; TEEND = 0.0 temperature during run DOS related values: ISMEAR = 0 ; SIGMA = 0.05 gaussian smear Electronic relaxation 2 (details) Write flags LWAVE = F write WAVECAR LCHARG = F write CHGCAR产生POSCAR和计算晶格常数的工作可以用以下的PBS script来完成。 #!/bin/tcsh #PBS -S /bin/sh #PBS -l nodes=4:athlon:ppn=2 #PBS -lcput=384:00:00 #PBS -m ae #PBS -o output #PBS -e error.log# set parameter set EXEC = vasp set SRC = /usr/common/executable# change working directory cd $PBS_O_WORKDIR# copy fresh executable from depository cp -f $SRC/$EXEC .# execute mpi program foreach a (3.3 3.4 3.5 3.6 3.7) echo a= $acat POSCAR SUMMARYend # remove executable rm -f $EXEC如果不用不需要用PBS script，则更加简单，如下即可。将其命名为lattice。 #!/bin/tcsh foreach a (3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2) echo a= $acat POSCAR SUMMARYend用chmod +x lattice，将其改为可执行文件。然后在命令行里键入./lattice 即可。 以下是用USPP-LDA运行完后的SUMMARY文件。每个计算用时13秒。（在USPP中Pd的截断能量是198.955） 3.5 1 F= -.52384500E+01 E0= -.52371846E+01 d E =-.253072E-02 3.6 1F= -.58695670E+01 E0= -.58683951E+01 d E =-.234381E-02 3.7 1 F=-.62322232E+01 E0= -.62311104E+01 d E =-.222547E-02 3.8 1 F=-.63932936E+01 E0= -.63921078E+01 d E =-.237151E-02 3.9 1 F=-.64072233E+01 E0= -.64058584E+01 d E =-.272979E-02 4.0 1 F=-.63162916E+01 E0= -.63147061E+01 d E =-.317085E-02 4.1 1 F=-.61523489E+01 E0= -.61504748E+01 d E =-.374817E-02 4.2 1 F=-.59418370E+01 E0= -.59396594E+01 d E =-.435530E-02用抛物线拟和得到的晶格常数为 ,固体中每个原子的能量是 。 以下是采用PAW-LDA势运行完以后的SUMMARY文件。每个计算用时20秒。所以相对来说PAW势所需要的时间多一些，这是因为PAW势的energy cutoff相对比较高（在PAW中Pd的截断能量是250.832）。 3.5 1 F= -.52393107E+01 E0= -.52377274E+01 d E =-.316665E-02 3.6 1F= -.58814938E+01 E0= -.58798653E+01 d E =-.325695E-02 3.7 1 F=-.62451262E+01 E0= -.62437004E+01 d E =-.285149E-02 3.8 1 F=-.64049388E+01 E0= -.64036223E+01 d E =-.263317E-02 3.9 1 F=-.64158100E+01 E0= -.64143798E+01 d E =-.286044E-02 4.0 1 F=-.63210060E+01 E0= -.63194198E+01 d E =-.317251E-02 4.1 1 F=-.61536329E+01 E0= -.61518107E+01 d E =-.364433E-02 4.2 1 F=-.59385695E+01 E0= -.59364165E+01 d E =-.430601E-02用抛物线拟和得到的晶格常数为 ,固体中每个原子的能量 可见，PAW-LDA和USPP-LDA给出的晶格常数都和实验吻合的非常好，两者之间的差别也很小。在以下所有的计算中，如果没有特殊声明，我们都默认采用PAW-LDA的势。结合能(cohesive energy)的定义如下： (104)单个Pd原子的能量为-1.426eV，所以我们得到Pd每个原子（相对于spin non-polorize的原子）的结合能为4.998eV。如果考虑Pd原子的spin-polarize的修正1.46eV，则结合能为6.458eV。用VASP计算表面能 做表面计算时，第一步我们需要测试K点的收敛性。通常，在垂直表面方向用1个K点就可以了，在平行表面方向，可以用和体材料类似的K点密度。 其次，我们要测试真空厚度(vacuum thickness)的收敛性。我们构造完一个slab后，将真空厚度逐渐从 增加到 ，体系的总能量改变不超过10meV的时候，可以初步认为真空厚度达到标准。以下是一个3层的(fcc) Pd slab的能量随着真空厚度的变化。其INCAR文件如下： SYSTEM = undeformed fcc Pd (111) surface calculation Startparameter for this run: PREC = Accurate ISTART = 0 job : 0-new 1-cont 2-samecut ICHARG = 2 charge: 1-file 2-atom 10-const ISPIN = 1 spin polarized calculation? Electronic Relaxation 1 NELM = 90; NELMIN= 8; # of ELM steps EDIFF = 0.1E-03 stopping-criterion for ELM LREAL = .FALSE. real-space projection NBANDS = 40 Ionic relaxation EDIFFG = 0.1E-2 stopping-criterion for IOM NSW = 0 number of steps for IOM IBRION = 2 ionic relax: 0-MD 1-quasi-New 2-CG ISIF = 2 stress and relaxation POTIM = 0.10 time-step for ionic-motion TEIN = 0.0 initial temperature TEBEG = 0.0; TEEND = 0.0 temperature during run DOS related values: ISMEAR = 1 ; SIGMA = 0.20 broadening in eV -4-tet -1-fermi 0-gaus Electronic relaxation 2 (details) Write flags LWAVE = F write WAVECAR LCHARG = F write CHGCAR LVTOT = .TRUE.其中因为Pd是金属，ISMEAR设置为method of Methfessel-Paxton。我们在最后的计算结果中必须保证entropy T*S这一项在OUTCAR中可以忽略不计（ ）。 POSCAR文件如下： Pd surface Calculation 3.875000000000000 0.7071067800000000 0.0000000000000000 0.0000000000000000 -0.3535533900000000 0.6123724000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 5.1961520000000000 4Selective Dynamics Direct 0.0000000000 0.000000000 0.0000000000 F F F 0.3333333333 0.666666667 0.1111111111 F F F 0.6666666667 0.333333333 0.2222222222 F F F 0.0000000000 0.000000000 0.3333333333 F F F如果对Direct的指定方法不熟的话，也可以用如下的POSCAR,和上面的完全等价。 Pd surface Calculation 1.0000 2.740038777 0.000000000 0.000000 -1.370019389 2.372943188 0.000000 0.000000000 0.000000000 20.135089 4Selective Dynamics Cartesian 0.0000000000 0.000000000 0.0000000000 F F F 0.0000000000 1.581962035 2.2372321073 F F F 1.3700193841 0.790981017 4.4744642247 F F F 0.0000000000 0.000000000 6.7116963320 F F FKPOINTS文件如下： fcc Pd K points 0 Monkhorst Pack 11 11 1 0 0 0下表列出了采用上面的参数设置，当真空层厚度从 增大到 时总能及功函数的变化： Table 4.1: 总能及功函数*随真空厚度的变化K points sampling : K points sampling : 真空厚度( 总能功函数(eV)总能功函数(eV)4-24.336866655.61006-24.213338015.62868-24.212872835.752510-24.21411465.724912-24.21277665.771714-24.21344785.856616-24.214554335.855618-24.212650955.870030-24.213636365.905550-24.213440035.6518*Pd(111)面功函数实验值为5.6。但是由于我们只用了4层原子，并且没有弛豫表面原子，所以并不能和实验直接对照。 我们可以看到，真空厚度大约为 时，体系的总能量就已经收敛。而如果要保证功函数的收敛，则真空厚度要加大到 左右。 首先我们要弛豫表面原子。弛豫的时候可以在INCAR中设置以下的参数。 POTIM=0.5 NSW=25 IBRION=2 EDIFFG=-0.03 MAXMIX=40另外，为了得到正确的结果，我们还需测试表面计算，特别是表面能是否收敛。我们必须保证实际计算中用到的slab模型足够厚，slab内部原子具有体材料的性质。如何判断slab内部原子具有体材料的性质呢？一个重要的标准就是当slab的层数N增加时，表面能的变化很小。 表面能的定义为 (105)A指的是每个表面上的原子数，2是因为我们有两个表面.表面能表示原子形成表面是所需要的能量，所以表面能越小的表面越稳定。 在slab计算中，一个很常用的用来计算表面能的公式是Boettger Equation(PRB 49:23) (106)这里， 指的是N层的slab的总能量。前面一个2指的是两个表面。后面两个2指的是stacking period。在逐渐增加slab层数的时候，还要注意同时保持超晶胞大小的一致。正如VASP手册上说的 It is almost impossible to compare two calculations which differ in the number of k-points and in the size of the supercell. 我们从3层开始一层层增加Pd的层数，以研究需要几层Pd原子才能达到收敛。在计算过程中，我们保持超晶胞大小的一致。 Table 4.2: 表面能收敛测试*K points sampling : K points sampling : N3-17.7636254-24.233309-6.4696841.64545-30.674313-6.4410041.53076-37.100998-6.4266851.459151.90377-43.560897-6.4598991.658451.89838-50.000495-6.4395981.516351.89879*在上述计算中，超晶胞大小不变。表面原子并未弛豫。包含弛豫后，表面能将降低。 采用VASP如何计算晶体的弹性常数弹性常数的概念#!RavindranP:Denftc:1998!#,#!Grimvall:tp:1999!# 弹性常数描述了晶体对外加应变 的响应的刚度。在应变很小的情况下，体系的内能与应变的大小存在二次线性关系（胡克定律），弹性常数 就是描述这种二次线性关系，即二次线性项的系数。采用Voigt标记： , , , , 和 。 应变张量 定义为： (107)应力张量 定义为： (108)二阶绝热弹性常数为： (109)在应变较小的情况下，应变后体系的总能 按应变张量 进可按泰勒级数展开为： (110)其中 是应变前体系的总能， 是应变前原胞的体积。另外，应变后基矢 与应变前的基矢 之间的关系为： (111)其中 为单位矩阵。 因此选取特定的应变 ，计算出在一组不同幅度时应变前后体系总能的变化( )，再根据总能的变化应变幅度对应的一组数据点，进行二次函数拟合得到二次项系数。即可得到晶体的某个弹性常数或弹性常数的组合。对不同的晶系的晶体，因为对称性的关系，它独立的弹性常数是确定的。比如对六角晶系的晶体，它独立的弹性常数为： , , , 和 。 六角晶体的弹性应变能与应变关系 对六角晶系的晶体，其原胞基矢可以取为： (112)其中 和 是晶体的晶格常数。 可以施加这样的应变 来计算 #!WangSQ:Abiec:2003!#： (113)施加应变 来得到 ： (114)对于 ，施加应变 来得到： (115)对于 ，施加应变 来得到： (116)最后，施加应变 得到 、 、 和 的组合： (117)由此可见，通过施加五个特定的应变，选取一系列幅度 的应变，得到 数据点，再分别按上面五个关系式对相应的 进行拟合得到二次项系数，最后联立方程得到求出六角晶系晶体的独立弹性常数。 具体计算步骤 在计算时，有几点要特别注意的：a)，原胞内原子在应变是否驰豫；b)，k点网格大小是否足够，因为在应变后，原胞的对称性会发生变化，即使同样的k点网格，在简约布里渊区产生的k点数目是不同的。因此，k点网格大小要取得足够，以保证弹性常数计算的精确性；c)，应变幅度 要取得适中，如果太小的话，得到的应变能（应变前后体系的变化）很小，在计算弹性常数时，会引起计算误差。 下面以计算六角AlN（纤锌矿结构）为例，在计算过程考虑了应变后原子的位置需要驰豫，具体步骤如下： 先对六角AlN体材料的晶格参数（晶格常数和原子位置）优化好，这样得到未应变时的POSCAR，并把它拷贝成下面defvector.f需要用到的一个的输入文件OLDPOS，并对OLDPOS做一个处理。因为OLDPOS在格式上有特殊要求： a 、在OLDPOS的第一行，在title的字符串之后，至少空一格再加上OLDPOS中原子的种类数目，比如AlN中有两类原子，title写为AlN，那么就OLDPOS的第一行就是AlN 2。 b 、OLDPOS的格式与POSCAR的类似，但是它最好是分数坐标来写出原子的位置。 以六角AlN为例，这个OLDPOS的内容如下： towidth 0.8pt height 0cm depth 0.25cmwidth 1.025height 0cm depth 0.8ptwidth 0.8pt height 0cm depth 0.25cm 1.0 AlN 2 3.11553 1.000000 0.000000 0.000000 -0.500000 0.866025 0.000000 0.000000 0.000000 1.605000 2 2Direct 0.00000000 0.00000000 0.00000000 0.333333333 0.666666667 0.50000000 0.00000000 0.00000000 0.381483673 0.333333333 0.666666667 0.881483673towidth 0.8pt height 0.25cm depth 0cmwidth 1.025height 0.8pt depth 0cmwidth 0.8pt height 0.25cm depth 0cm 对特定的应变，在下面的defvector.f中Define the strain部分，把特定的应变通过给strain(i)矩阵赋值。比如对上面提到的 用来计算 ，那么就对defvect.f中的Define the strain改写成如下的形式： towidth 0.8pt height 0cm depth 0.25cmwidth 1.025height 0cm depth 0.8ptwidth 0.8pt height 0cm depth 0.25cm 1.0 C% Define the strain % strain(1)=delta strain(2)=delta strain(3)=0.0 strain(4)=0.0 strain(5)=0.0 strain(6)=0.0towidth 0.8pt height 0.25cm depth 0cmwidth 1.025height 0.8pt depth 0cmwidth 0.8pt height 0.25cm depth 0cm 其中整个defvector.f是用来得到某个应变后，新的POSCAR。应变的类型按Define the strain的部分来定义，而 应变的幅度需在程序编译后，运行编译得到的模块时输入。defvect.f的内容如下： towidth 0.8pt height 0cm depth 0.25cmwidth 1.025height 0cm depth 0.8ptwidth 0.8pt height 0cm depth 0.25cm 1.0 C%C this simple program to get the primitive vectors after C$delta$ strain, in order to calculate the independent Celastic constants of solids. C usage: C! Please firstprepare the undeformed POSCAR in OLDPOS C defvector.x Ctype defvector.x create new POSCAR in file fort.3C% program defvector real*8 privect,strvect,delta,strten,strain,pos, alat dimension privect(3,3),strvect(3,3),strten(3,3),strain(6) dimension pos(50,3) character*10 bravlat, title, direct integer i,j,k,ntype, natomi, nn dimension natomi(10)C% Read the undeformed primitive vector and atomic postion % open(7,file=OLDPOS)C% In first line of OLDPOS, please add the numberC% of the type of atoms after the title read(7,*) title, ntype read(7,*) alat do i=1,3 read(7,*) (privect(i,j),j=1,3) write(*,*) (privect(i,j),j=1,3) enddo read(7,*) (natomi(i),i=1,ntype) nn=0 do i =1, ntype nn=nn+natomi(i) enddo read(7,*) direct do i=1, nn read(7,*) (pos(i,j),j=1,3) enddoC% Read the amti of strain % read(*,*) deltaC% Define the strain % strain(1)=delta strain(2)=0.0 strain(3)=0.0 strain(4)=0.0 strain(5)=0.0 strain(6)=0.0C% Define the strain tensor % strten(1,1)=strain(1)+1.0 strten(1,2)=0.5*strain(6) strten(1,3)=0.5*strain(5) strten(2,1)=0.5*strain(6) strten(2,2)=strain(2)+1.0 strten(2,3)=0.5*strain(4) strten(3,1)=0.5*strain(5) strten(3,2)=0.5*strain(4) strten(3,3)=strain(3)+1.0C% Transform the primitive vector to the new vector under strain%C strvect(i,j)=privect(i,j)*(I+strten(i,j) do k=1,3 do i=1,3 strvect(i,k)=0.0 do j=1,3 strvect(i,k)=strvect(i,k)+privect(i,j)*strten(j,k) enddo enddo enddoC% Write the new vector under strain% do i=1,3 write(*,100)(strvect(i,j),j=1,3) enddo 100 format(3f20.15)C% Create the POSCAR for total energy calculation %5 write(3,(A10) title write(3,(f15.10) alat do i=1,3 write(3,100)(strvect(i,j),j=1,3) enddo write(3,(10I4) (natomi(i), i=1,ntype) write(3,(A6) Direct do i=1, nn write(3,100) (pos(i,j),j=1,3) enddoC% endtowidth 0.8pt height 0.25cm depth 0cmwidth 1.025height 0.8pt depth 0cmwidth 0.8pt height 0.25cm depth 0cm 在defvector.f中，原子种类的数目由变量ntype来定义的，最大为10，如果原子的种类数目太大，需自己手动调大数组natomi(10)以及输出时 write(3,(10I4) (natomi(i), i=1,ntype)的格式10I4。 在defvector.f中设置好了特定的应变后，就可以编译defvector.f（使用g77 -o defector.x defector.f)得到模块defvector.x。 准备好VASP计算的输入文件KPOINTS和POTCAR。以及进行原子位置驰豫计算的INCAR.relax，它的内容如下： towidth 0.8pt height 0cm depth 0.25cmwidth 1.02