一、用VASP计算态密度

使用的软件：VASP

辅助分析计算的小程序：read_dos.f90

操作流程及主要步骤：

主要分两步：一、静态自洽计算；二、非自洽计算（以fcc结构的Al为例）

第一步静态自洽计算（得到自洽的电荷密度CHG、CHGCAR和E-fermi，提供给下一步非自洽计算用）

• 准备INCAR：即定义PREC，EDIFF，ENCUT，ISTART= 0（默认），IBRION = 2，ISMEAR(对于半导体和绝缘体：ISMEAR= 0，SIGMA = 0.05；ISMEAR= 0 表示采用Gaussian smearing 方法； 对于金属：ISMEAR= 1，SIGMA = 0.2;ISMEAR= 1 表示采用1阶smearing方法)

                   SYSTEM=Al-fcc

                   PREC=Accurate

                   EDIFF= 1e-5

                   ENCUT= 250.0

                   ISTART= 0 ；IBRION = 2  #【默认NSW=0，此处定义IBRION=2无意义】

           ISMEAR= 1 ; SIGMA = 0.2

• 准备好POSCAR文件，或以优化的晶格参数作为基础，把结构优化产生的CONTCAR拷贝成POSCAR

                   Al

                     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

• POTCAR直接从赝势库中得到

• KPOINTS

                Automatic generation

                0

               MohkorstPack

                  9    9    9  

                  0.0  0.0 0.0

• 提交运行

• 记下OUTCAR中E-fermi的数值，下一步处理结果时会用到。

              grep“E-fermi”OUTCAR

第二步非自洽计算

• 准备INCAR：即定义ENCUT，PREC，ISTART= 1（存在WAVECAR文件,所以取1），ICHARG = 11（表示从CHGCAR中读入电荷分布，并且在计算中保持不变），ISMEAR= -5 ( -5表示带Bloch修正的四面体方法；采用 tetrahedron method with Bl¨ochl corrections计算DOS更准确一些，在PWscf中也是如此) ，RWIGS (或LORBIT = 10，这时可不设RWIGS)
  

SYSTEM = Al-fcc

PREC = Accurate

ENCUT = 250              #【截断能不能大于上一步非自洽计算？最好与上一步计算保持一致】        

ISTART = 1；ICHARG = 11

ISMEAR = -5

     LORBIT = 10

# 【对于PAW势,可设置LORBIT = 10，此时可不用设置RWIGS参数；具体解释见附】

#【一般上述参数用来做非自洽计算、产生DOS已经足够；如NSW表示离子运动步数，默认为0(此时，IBRION默认为-1，表示 原子位置不移动)，即做静态计算。自洽还是非自洽，如何判断？】

• 准备好KPOINTS文件，增加k点网格(增加k点可加大积分计算精度，因为是非自洽计算，k点数虽增大但仍可以很快)

      Automatic generation

      0

     MohkorstPack

         19   19  19      

          0.0  0.0 0.0

• POSCAR、POTCAR同上一步自洽计算。

• 将上一步自洽计算得到的CHG、CHGCAR拷贝至同一目录下，提交运行。【ICHARG=11，ISMEAR=-5】

• 计算完后，得到包含了态密度值的DOSCAR文件，采用read_dos.f90小程序对态密度文件DOSCAR进行处理，得到总态密度dos_total.dat，及分波态密度*_dos.dat。

   运行read_dos.f90时会提示：please input the Ef，此时需要把第二步自洽计算得到的E-fermi的值给出。接着会提示：please input the number of atoms in the POSCAR file，此时需要把POSCAR中的原子数给出。【操作./read_dos.x  DOSCAR】

   输出文件的能量值是以费米能级作为能量参考零点。dos_total.dat的第一列数据是能量值，单位为eV；第二列数据是总态密度的值，单位State/eV.unit cell；第三列数据是总态密度的积分值，也就是电子数，单位为electrons (用origin等软件绘图时把dos_total.dat的最后一列删除即可) 。    

     *_dos.dat是相应原子【unit cell中含有的每个原子都有各自对应的s、p、d分波态密度值】的分波态密度值，其中的第一列数据是能量值，单位为eV；第二、三、四列数据分别对应于s、p、d态的分波态密度值，单位为State/eV.atom。

#【如果将s、p、d态密度按每种元素求和，则态密度和的单位是State/eV.atom?如何将单位转换为State/eV.el.】

#【为什么自洽完了，要做非自洽计算的英文解释见后面附中的英文解释。提问：非自洽产生的DOSCAR 文件与上一步自洽计算产生的DOSCAR之间有什么区别？】

二、用VASP计算能带

使用的软件：VASP

辅助分析计算的小程序：read_band.f、get_kpoints.sh

操作流程及主要步骤：

计算材料的能带结构即色散曲线E（k），主要分两步：一、静态自洽计算(fcc)；二、非自洽计算（以fcc结构的Al为例）

第一步静态自洽计算（同上）【因此，可单独建一个scf的文件夹，用以产生并存储电荷密度CHG、CHGCAR和E-fermi】

第二步非自洽计算（在固定电子密度的情况下，得到选取K点的能量本征值。）

• 准备INCAR：即定义ENCUT，PREC，ISTART= 1，ICHARG = 11（表示从CHGCAR中读入电荷分布，并且在计算中保持不变），并增加NBANDS（默认值为NELECT/2+NIONS/2，NELECT和NIONS分别为电子数和离子数，可以上一步自洽计算产生的OUTCAR文件中找到这两个参数，如grep "NIONS" scf/OUTCAR，和 grep "NELECT" scf/OUTCAR。)【NBANDS的设置，详见后面的英文注释】

SYSTEM = Al-fcc

PREC = Accurate

EDIFF = 1e-5

ENCUT = 250

IBRION= 2

ISTART = 1；ICHARG = 11

ISMEAR= 1；SIGMA =0.2

NBANDS = 12

• 准备KPOINTS文件，使用Line-mode，给出高对称性k点之间的分割点数。【分割越密，则路径积分越准确，VASP 5.x格式有变吗？】

      k-points along high symmetry lines    !注释行，无特别的意义

          20                                                  ! intersections，沿G-X特殊点之间产生10个k点

     Line-mode                                       !  程序自动产生特殊k点间的k点

     Rec                                                 ! 各k点相对于倒格子基失来写的  Al的元胞的高对称性点

       0.5      0.0         0.5   ! X

       0.0      0.0         0.0   ! G

      0.0      0.0         0.0   ! G

      0.5      0.5         0.5   ! L

     0.5      0.5         0.5   ! L

     0.5      0.25        0.75  ! W

     0.5      0.25        0.75  ! W

     0.375   0.375       0.75  ! K

     0.375   0.375       0.75  ! K

     0.0      0.0         0.0    ! G

#【对高对称性点的理解及注释见后附注部分。】

说明：通过指定Line-mode，VASP会自动在起点和终点之间插入指定的K点数，比如上面的文件就是指定VASP计算沿着X到Gamma点再到L、W、K回到Gamma点，每个方向上各取20个K点。

• POSCAR、POTCAR同上一步自洽计算。

• 将上一步自洽计算得到的CHG、CHGCAR拷贝至同一目录下，提交运行，计算完后可以从OUTCAR文件或者EIGENVAL文件里得到需要的每个K点的能级信息【与运用WAVECAR得到的能带结构一样，CHGCAR是WAVECAR的积分。是的，波函数代表一切】  

• 运行小程序get_kpoints.sh，得到kpoints.dat（即高对称点的位置；）。编译程序read_band.f，运行，系统提示：

     please input thefermi energy:E-fermi from scf OUTCAR

      将上一步得到的E-fermi值输入到屏幕，则输出bands.dat，并产生spoint.dat文件 (该文件给出了作图时高对称性点的横坐标位置)，因此，可用origin等绘图软件作图了。【read_band.f程序的解释见附】

lpf的经验和心得：

  对于第一步自洽计算：由于在能带计算时k点是一些在倒空间高对称线上的点，不能进行自洽计算【解释得好啊】，因此在进行能带计算前必须加一次自洽计算以得到精确的电荷密度值。

   自洽计算得到的电荷密度文件CHGCAR是能带和DOS计算需要的输入文件。

   对于非自洽计算：计算能带时ICHARG= 11，金属用ISMEAR=1；半导体或绝缘体，用ISMEAR=0 。

   计算态密度：ICHARG = 11，ISMEAR=-5

     【WAVECAR文件在计算DOS和BAND过程中有用么？是万能的，但要开计算band和dos的源程序调用谁了】

   用origin绘图技巧，用Vi打开spoint.dat文件，得到作图时高对称点的横坐标位置，将横坐标端点分别设为第一个、最后一个高对称点对应的坐标；画竖线，双击，修改Objectproperties---Cordinates:Units选为Scale，设置开始x和终点x均为spoint.dat文件中的高对称点数值；y值即设为纵坐标范围。【实际上图形中的横坐标范围已经定死了，无法修改；其实也恰好与spoint.dat中点的范围一致。如何断定spoint.dat文件中数值具体对应哪个高对称点符号？】

1：赝势

     赝势文件夹各个目录对应起来分别是pot ==> PP, LDA ;                                 pot_GGA ==> PP, GGA ;        potpaw ==> PAW, LDA ;

                                                             potpaw_GGA ==> PAW, GGA, PW91 ;   potpaw_PBE ==>PAW , GGA, PBE。

      选择某个目录进去，我们还会发现对应每种元素往往还会有多种赝势存在。这是因为根据对截断能量的选取不同还可以分为Ga,Ga_s,Ga_h，或者根据半芯态的不同还可以分为Ga,Ga_sv,Ga_pv的不同。

2：non-scf

      ISTART=1;ICHARG=11时，PREC必须与做静态计算的设置一致（尤其是截断能），否则会报错

      WARNING:dimensions on CHGCAR file are different

      ERROR: chargedensity could not be read from file CHGCAR for ICHARG>10

3：POSCAR

     POSCAR中的那个缩放系数对于正格矢起作用如果原子坐标是用Car格式，它也起作用。如果原子坐标是Direct格式，则它不起做作用。【yexq评：缩放系数即对应PWscf中的alat，只影响晶格基矢量的大小，个人感觉不会影响原子的crystal位置】。

4：INCAR

    NSW是优化原子位置和晶胞个开启项，默认=0。如果没有设置NSW，即使ISIF=3，依然只做静态计算。

5：KPOINTS

     采用Line-mode时，”0 0 0 ！gama“的最后一个0和！之间需要有空格。【yexq评：真是经验总结，很仔细】

附录-yexq总结：

1.    为什么自洽计算完了，还要做非自洽计算，请见下面的英文解释

     Calculating a DOS can be done in two ways:

     The simple one is toperform a static (NSW=0, IBRION=-1) selfconsistent calculationand to take the DOSCAR file from this calculation.    

     However, the simple approach discussed above is not applicable in all cases: A high quality DOS requires usually very fine kmeshes.  You should think at least in orders of 16x16x16 meshes for small cells and even for large cells you might need something like 6x6x6- or 8x8x8-meshes.

      The usual way,to do DOS or band-structure calculations in this case is the following: the charge density and the effective potential converge rapidly with increasing number of k-points.

     So, as a first step one generates a high quality charge density using a few k-points in a static selfconsistent run. The next step is to perform a non-selfconsistent calculation using the CHGCAR file from this selfconsistent run (i.e. ICHARG is set to 11, see section 6.14) .

     Mind, this is the only way to calculate  the band structure, because for band-structure calculationsthe supplied k-points form usually no regular three-dimensional grid and therefore a selfconsistent calculation gives pure nonsense ! 【但是还是需要前面的自洽计算】

      For ICHARG=11, all k-points can be treated independently, there is no coupling between them, because the charge density and the potential are kept fixed. Therefore there is also no need to treat the k-points within one single run simultaneously.

   【可新建1个pdos文件夹，非自洽计算可安排在该文件夹中进行】

2. NBANDS的设置问题

   可见NELECT是所采用赝势中的价电子数吗？例如 对TiH2：取14/2 +  3/2 =  8.5，PdH： 取17/2+2/2=9.5

3. 对高对称性点选取的理解

     采用与自洽计算所用的相同大小的胞(元胞或单胞)，来选取高对称性点；

   高对称性点并不一定形成闭合的曲线；

   应对结构施加正确的symmetry后，选高对称性点。这也说明，能带计算不是在uniform的格子上进行，而是沿高对称性点划定的路径积分。所以，如果结构没施加对称性，则选取的高对称性点将与施加了对称性操作的“相同”结构不一样；同理按元胞找的对称性点与按单胞找的对称性点将不一样。

   对于具有相同对称性的结构，其高对称性点都有一致的约定的取法，即与具体物质无关。例如同为fcc结构的不同物质，其高对称性点将是一样的。因此可以通过下述网站输入结构的空间群号寻找高对称性点：

  http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-kv-list，或采用Material studio 得到高对称性点。

4.能带计算完成后，绘制能带图时，如何确定图中高对称性点的位置。    

   spoint.dat文件中数值按大小排列，其对应的高对称性点的顺序可能与计算前KPOINTS中输入的高对称点的顺序一致，实际上高对称性点的crystal坐标(即KPOINTS中的三个数值)可在EIGENVAL文件中找到(或OUTCAR中找到)。

    其实，read_band.f程序很简单，即将路径上的各k点crystal坐标求均方根，然后逐点累加，并记路每一k点对应的engenvalue；当遇到特殊k点(高对称性点)时，同时将该点坐标输出到spoint.dat文件中。  

     该程序跟据KPOINTS文件中特殊k点之间的取样间隔来判断某点是否是特殊k点，如是则同时输出其坐标到到spoint.dat文件中。因此，当特殊k点之间插入的一般k点数变化时，则应修改源文件中的if(mod((i-1),20).eq.0 .or.i.eq.nks语句，例如将20修改为所插入的一般k点数或实际分割的k点数。因此，该程序缺乏通用性，然而程序虽简单，但却高度体现了作者对能带结构的理解。

    实际高对称性点的位置可以直接求均方根确定吗？(不能，没有加和，还应该加上路径上其他点坐标产生的加和)

2.部分参数的英文注释

（1）ISMEAR = -5 | -4 | -3 | -2 | 0 | N

SIGMA = width of the smearing in eV

Default

ISMEAR = 1

SIGMA = 0.2

ISMEAR determines how the partial occupancies fnk are set for each wavefunction. For the finite temperature LDA SIGMA determines the width of the smearing in eV.

ISMEAR:

-5 tetrahedron method with Bl¨ochl corrections

    For the calculation of the total energy in bulk materials we recommend the tetrahedron method with Bl¨ochl corrections (ISMEAR=-5). This method also  gives a good account for the electronic density of states (DOS).

       The only drawback is that the methods is not variational with respect to the partial occupancies. Therefore the calculated forces and the stress tensor can be wrong by up to 5 to 10 % for metals.

       For the calculation of phonon frequencies based on forces we recommend the method of Methfessel-Paxton (ISMEAR>0). For semiconductors and insulators the forces are correct, because partial occupancies do not vary and are zero or one.

     For the calculations of the DOS and very accurate total energy calculations (no relaxation in metals) use the tetrahedron method (ISMEAR=-5).

      For semiconductors or insulators use the tetrahedron method (ISMEAR=-5), if the cell is too large (or if you use only a  single or two k-points) use ISMEAR=0 in combination with a small SIGMA=0.05.

     For relaxations in metals always use ISMEAR=1 or ISMEAR=2 and an appropriated SIGMA value (the entropy term should be less than 1 meV per atom). Mind: Avoid to use ISMEAR>0 for semiconductors and insulators, since it might cause problems. For metals a sensible value is usually SIGMA= 0.2 (which is the default).

（2）RWIGS or LORBIT

    If RWIGS or LORBIT (Wigner Seitz radii, see section 6.326.33) is set in the INCAR file, a lm- and site-projected DOS is calculated and also written to the file DOSCAR. One set of data is written for each ion, each set of data holds NDOS lines with the following data

             energy s-DOS p-DOS d-DOS

or

              energy s-DOS(up) s-DOS(down) p-DOS(up) p-DOS(dwn) d-DOS(up) d-DOS(dwn)

for the non spin-polarized and spin polarized case respectively. As before the written densities are understood as the difference of the integrated DOS between two pins.

(3) DOSCAR file

    Mind: For relaxations, the DOSCAR is usually useless. If you want to get an accurate DOS for the final configuration, first copy CONTCAR to POSCAR and continue with one static (ISTART=1; NSW=0) calculation.【静态自洽计算？】

（4）LORBIT

            The default for LORBIT is .FALSE. (respectively 0).

           This flag determines, together with an appropriate RWIGS (see section 6.32), whether the PROCAR or PROOUT files (see section 5.21) are written.

           The file PROCAR contains the spd- and site projected wave function character of each band.

            The wave function character is calculated, either by projecting the wavefunctions onto spherical harmonics that are non zero within spheres of a radius RWIGS around each ion (LORIT=1, 2), or using a quick projection scheme relying that works only for the PAW method (LORBIT=10,11,12, see below). If the LORBIT flag is not equal zero, the site and l-projected density of states is also calculated.

           If the projector augmented wave method is used, LORBIT can also be set to 10, 11 or 12. This alternative setting selects a quick method for the determination of the spd- and site projected wave function character and does not require the specification of a Wigner-Seitz radius in the INCAR file (the RWIGS line is neglected in this case). The method works only for PAW  POTCAR files and not for ultrasoft or norm conserving pseudopotentials.  即对每个能带的波函数进行s  pd和site分解或投影；LORBIT = 10，不设RWIGS，针对PAW势的计算。

   LORBIT = 10表示对波函数按原子site进行s 、p、d轨道分解或投影；并输出投影波函数到PROCAR or PROOUT 文件中，详见上面英文注释部分：This alternative setting selects a quick method for the determination of the spd- and site projected wave function character。什么样的quick method?