Martin's profileMartin's FactoryPhotosBlogLists Tools Help

Blog


    September 26

    赝势方法;模型势方法;pseudo-potential method

     分子式:
    CAS号:

    性质:又称模型势方法。即价电子从头计算法。其基本思想是,所研究体系的哈密顿算符仅显含价电子部分,而将原子内层的全部电子连同原子核构成的核实对外层价电子的作用用适当的模型势函数(即赝势pseudo potential)表示,同时引入表示投影算符的势以便将价电子波函数与内层电子波函数分离开来,然后对价电子进行变分并用自洽迭代处理,计算出价电子轨道波函数和能级值等。所引入的两个势均通过全电子的原子从头计算确定,随后用于分子计算。常见的赝势有两类,一类是基于固体物理理论的菲利普斯-克兰曼(Phillips-Kleinman)赝势原理发展起来的,一类是从早期的价电子概念出发,由赫律纳加(Huzinaga)等提出的模型。赝势方法主要用于含有重元素的化合物或原子簇的研究。

    VASP赝势文件POTCAR详解

     
    POTCAR 赝势文件

      可以理解为分子力学模拟中的力场文件 但包括的信息更多
      VASP4.6将各元素优化的INCAR里的参数也包括在这里了,作为支持PREC的缺省选择
      通常各元素的POTCAR已经包括在软件包里了
      我们只需要按照POSCAR里的顺序,将各元素的POTCAR按顺序连接起来就可以了
      如以下命令:
      cat file1 file2 file3 > POTCAR

      软件包自带的绝大多数赝势是超软赝势(US-PP)了,但不少元素有两个版本,如何
      选取呢?
      一个简单的办法是看后缀
      标准的没有后缀 _h 硬一点 _s 软一点
      _pv,_sv,_d 就是说semi-core的p,s或者d也当做价态处理了
      如果是数字的话,表示的可能是不同的半径截距

      也可以参考各版本同目录下的V_RHFIN file ,PSCTR file
      这两个文件告知该版本的赝势是如何生成的。比如:
      V_RHFIN file
      Sc: 6p d2 s1
       8 21. .002000 44.95590 125. .25E-05 .300 200FCA 12.00000
       .7 1.0 0
       1.0 .0 .5 -320.8847 2.0000
       2.0 .0 .5 -34.4217 2.0000
       2.0 1.0 1.5 -28.2366 6.0000
       3.0 .0 .5 -3.7944 2.0000
       3.0 1.0 1.5 -2.2591 6.0000
       3.0 2.0 2.5 -.1113 2.0000
       4.0 .0 .5 -.2699 1.0000
       4.0 3.0 2.5 -.1000 .0000
      第一行是注释行 给出基本的信息
      第二行是最重要的控制行
       8 21. .00 2000 44.95590 125. .25E-05 .300 200 F CA 12.00000
       J Z XION N AM H DELRVR PHI NC1 | CH QCOR
       |
       GREEN

      J - 轨道数 Z - 原子序数 XION - 离子化程度 一般设为0 N - 格点数
      AM - 原子质量 H - 决定格点间距 DELRVR - 自洽收敛标准
      PHI - 线性拟合参数 NC1 - 最大自洽循环次数 GREEN - 是否存在初始的势
      CH - 交换相关能(XC)类型
       Slater-XC
       HL Hedin Lundquist (1971)
       CA Ceperly and Alder parameterized by J.Perdew and Zunger
       WI Wigner interpolation
       PB Perdew -Becke
       PW Perdew -Wang 86
       LM Langreth-Mehl-Hu
       91 Perdew -Wang 91
      QCOR - 非价键电子数(core electrons)
      第三行开始是每个轨道的具体参数,依次为
      n l j(=l±1/2) 原子轨道能 占有率

      PSCTR file of LDA/H1.25
       TITEL = US H
       LULTRA = T use ultrasoft PP ?
       RWIGS = 0.57 nn distance ! Wigner-Seitz radius

       RCLOC = .65
       NE = 100
       LCOR = .TRUE.
       QCUT = -1
       RMAX = 3.0 ! core radius for proj-oper

       Description
       l E TYP RCUT TYP RCUT(cutoff radius)
       0 0 15 0.80 23 1.25
       0 0.5 15 0.80 23 1.25
       1 -0.2 15 0.80 23 1.25
      最重要的地方上面已经用颜色标出来啦:)
      说明一下,TYP是指赝势的类型,RCUT是半径截距,TYP可取的值如下:
      正则
      1 BHS
      2 TM
      3 VAN
      6 XNC
      7 RRKJ wave function possibly with node
      15 RRKJ wave function strictly no node
      非正则 +8

      最后一个问题是LDA or GGA。貌似没有定论目前。
      这个最好是两个一起做做看啦。或者看文献别人验证过哪个数据好。
      其实据说目前最好的是PAW(P.E.Blochl,Phys.Rev.B 50,17953(1994).,Phys.Rev.B 59,1758(1999).)。

    VASP新手入门十个简单问题

    1.TEBEG,TEEND的单位是K吗?
      答:是K为单位的。

      2.在运算中,直接COPY POSCAR文件到工作目录就可以吗,用不用再POSCAR或INCAR文件中设定参数来控制POSCAR的赝势的选取方法?
      答:不大明白你的问题。在一个工作目录准备好INCAR,POSCAR,POTCAR,KPOINTS这四个文件。先看POTCAR中LEXCH后面的值是什么,如果是CA,那么INCAR中就不用设置GGA了;如果是91,那么就在 INCAR 中输入GGA=91;如果是PE,那么就设置GGA=PE。

      3.模拟表面,一般SLAB选多厚?
      答:自己多测试不同的厚度后,来确定slab 的厚度

      4.在KPOINT文件中,如果选取MONKHORST-pack的K点产生办法,K点相对自动网格移动时怎么回事?
      答:你可以看看monkhorst,pack1976年在phys. rev. b上的文献,他们在文献中给出了对于sc,bcc和fcc三种格子的一种效率很高的k点选取标准。但是这种k-mesh的选取方法并不是将第一个点取在倒空间的原点处,而是有一个shift,比如对于fcc而言,这个shift时(1/8,1/8, 1/8)。但是这种mp方法对于正交格子适用,而对于三方或者六方的格子则并不好。因此这时建议采用\gamma点为中心的k-mesh(vasp用户手册)

      
      5.想做表面重构,VASP 能做吗,能用分子动力学模拟吗,怎么设置参数?
      答:做重构的计算,一般用分子动力学的方式来模拟。你可以参考文献上模拟表面体系重构的作法。大多是采用比较几个表面模式的能量以及相关的热力学量。
      
      6.不指定电子的总SPIN,会有影响吗,一般 什么体系才考虑自旋?
      答:你做的体系是否包含3d-过渡金属离子(V至Ni)或其他含f电子的原子,如果没有的话,可以不考虑自旋极化的计算。
      
      7.表面的原子或分子的运动,扩散还有重构,IBRION怎么设定?
      答:是扩散的化,一般用NEB的方法计算其扩散路径或扩散势垒。

      8.如果知道一个化合物的组成,想计算它的总能一般应怎么办?同时如果一直组成如何在常温下计算它的稳定结构,比如FeMgCl2?
      答:计算总能,自己准备好vasp的输入文件,运算后就直接得到。可能考虑它的结合能比较好。

      9.在计算一个系统的稳定结构时,SMASS参数如何设置?
      答:看你是计算这个稳定结构的什么性质:计算总能和DOS,用ISMEAR=-5,不用管SMASS。

      10.具体给介绍一下STM图像的模拟办法,如何计算特定的能带,还有计算完成后一般去那个文件查找对自己有用的数据?
      答:一般来说用parchg做个图就可以了,当然以上方法仅限于T-H近似。 

    如何用VASP计算晶格常数

     

    我们用Pd金属作为例子。
      Pd金属的实验上的晶格常数为3.89A。在这里,我们用VASP计算它的晶格常数。
      首先将Pd所对应的POTCAR文件拷贝到目录下。然后准备好INCAR和KPOINTS文件。POSCAR文件我们将通过一个tcsh的script来产生。

      KPOINTS文件可以如下:
      Monkhorst Pack
      0
      Monkhorst Pack
       11 11 11
       0 0 0

      INCAR文件可以如下:
       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 -l cput=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= $a"

      cat >POSCAR <
      cubic diamond
       $a
       0.5 0.5 0.0
       0.0 0.5 0.5
       0.5 0.0 0.5
       2
      direct
       0.0 0.0 0.0
       0.25 0.25 0.25
      !

      mpiexec -nostdin ./$EXEC

      set E=`tail -2 OSZICAR`
      echo $a $E >>SUMMARY

      end
      # 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= $a"

      cat >POSCAR <
      fcc lattice
       $a
       0.5 0.5 0.0
       0.0 0.5 0.5
       0.5 0.0 0.5
       1
      cartesian
       0.0 0.0 0.0
      !

      ./vasp

      set E=`tail -1 OSZICAR`
      echo $a $E >>SUMMARY

      end

      用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 1 F= -.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
      用抛物线拟和得到的晶格常数为$3.888\AA$,固体中每个原子的能量是$E_{bulk}=-6.4257$。

      
      以下是采用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 1 F= -.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
      用抛物线拟和得到的晶格常数为$3.875\AA$,固体中每个原子的能量E_bulk=-6.4185eV

      
      可见,PAW-LDA和USPP-LDA给出的晶格常数都和实验吻合的非常好,两者之间的差别也很小。在以下所有的计算中,如果没有特殊声明,我们都默认采用PAW-LDA的势。
      结合能(cohesive energy)的定义如下:
      -E_coh = [E_bulk-N*E_atom]/N
     
    所以我们要将固体中每个原子的能量减去单个Pd原子的能量,才能和实验的结合能相比较。对于过渡金属原子,计算单个原子的能量要特别注意。VASP的网页上给出了求结合能所需的单个原子能量的修正值(详见VASP手册Pseudopotentials supplied with the VASP package一章)。可以在上面查到,Pd每个原子LDA的修正值为1.46eV。所以我们得到LDA近似下Pd的结合能为4.998eV。此值和实验值比严重偏大,这是因为LDA
    通常成键过强的关系。如果我们改用GGA的赝势,可以得到和实验比较吻合的结果。  

    如何用VASP计算单个原子的能量和能级

     氢原子的能量为-13.6eV在这一节中,我们用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
       1
      cart
       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

      
      我们可以看到,电子的能级不为-13.6eV。

      
       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

      
      我们可以看到TOTEN-EATOM也不等于-13.6eV。

      在上面的计算中有个问题,就是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

      氢原子的能量约等于-12.92eV。可以看到在LDA中如果限制自旋,使能级大概提高了0.88eV。

      如果我们采用GGA的赝势,并且同样打开自旋限制(ISPIN=2),在此例子中,得到的结果将更加精确
       Free energy of the ion-electron system (eV)
       ---------------------------------------------------
       alpha Z PSCENC = 0.00621465
       Ewald energy TEWEN = -1.38027565
       -1/2 Hartree DENC = -6.91107031
       -V(xc)+E(xc) XCENC = 2.32601856
       PAW double counting = 0.70286470 -0.71256934
       entropy T*S EENTRO = 0.00434843
       eigenvalues EBANDS = -7.71257941
       atomic energy EATOM = 12.52153358
       ---------------------------------------------------
       free energy TOTEN = -1.15551480 eV

       energy without entropy = -1.15986323 energy(sigma->0) = -1.15696428

      用GGA算得的氢原子的能量约等于-13.67eV。

      但是如何理解所得到的能级,由于用到了赝势,本人并不很清楚如何解释能级意义,欢迎大家指教。

    如何用VASP计算单个原子的能量[权威版]

     这是VASP作者给出的解答,所以算是权威了.原文如下:

      The energy zero
      ..
      in VASP all energies are referred to the the reference state for which
      the potential was generated!
      this is in most cases not the real groundstate of the atom
      ..
      to determined the energy of the grounstate of the atom place the atom
      in a larger non cubic box to break initial symmetry (i.e.
      11Ax10Ax9A)
      use the G point only
      INCAR:
      ISPIN = 2 ! spin polarized
      ISMEAR = 0 ; SIGMA = 0.2 ! for small sigma conv. for TM is diff.
      MAGMOM = 2 ! initial magnetic moment
      one should use the energy value energy without entropy of the OUTCAR file
      since this coverges most rapidly to the correct energy for sigma 0
      ..
      
       Ecoh=E_total-n*Eatom

      大约翻译如下:
      在VASP中,所有的能量是相对于产生赝势的组态来说的,所以并不能将计算得到的能量当成真正的原子的基态能量。

      为了得到原子的基态能量,将原子放在一个非正方体的晶格内,比如说11Ax10Ax9A的晶格(这样做是为了消除简并)。
      KPOINTS只用G撒点
      大概的INCAR文件可以如下:
      ISPIN = 2 ! spin polarized
      ISMEAR = 0 ; SIGMA = 0.2 ! for small sigma conv. for TM is diff.
      MAGMOM = 2 ! initial magnetic moment

      计算结束后,应该用OUTCAR中不含entropy的那一项能量,因为它通常更快的收敛到sigma=0的正确值。最后计算结合能的公式应该为
      
       Ecoh=E_total-n*Eatom

    用VASP计算氢分子的能量

     氢气分子的解离能,也就是结合能,根据资料中给出的是约4.48eV。(G. Kresse & J. Hafner, Surface Sci. 459 (2000) 287)
      
      为此,首先要计算一个氢气分子的孤立能量,再减去两个孤立氢原子能量,将得到氢气分子的结合能。
      
      先计算单个原子能量,选取PAW_PBE文件夹下的H下面的POTCAR
      用到的其他输入文件如下:
      INCAR:
      SYSTEM = H atom in a box
      ISMEAR = 0 ! Gaussian smearing
      SIGMA = 0.01
      ENCUT = 350.0

      KPOINTS:
      Automatic mesh
       0
       Monkhorst Pack
       1 1 1
       0. 0. 0.
      
      POSCAR:
      H atom in a box
      1.0 ! universal scaling parameters
      7.0 0.0 0.0 ! lattice vector a(1)
      0.0 8.0 0.0 ! lattice vector a(2)
      0.0 0.0 9.0 ! lattice vector a(3)
      1 ! number of atoms
      cart ! positions in cartesian coordinates
      0 0 0

      用上述文件计算得到 TOTEN(H) = +0.000854eV (H atom),用PAW_GGA得到类似结果。
      
      然后计算氢分子能量,用类似的输入文件:
      INCAR:
      SYSTEM = H2 dimer in a box
      ISMEAR = 0 ! Gaussian smearing
      NSW = 5 ! 5 ionic steps
      IBRION = 2 ! use the conjugate gradient algorithm
      ENCUT = 350.0
      POTIM = 0.1

      KPOINTS:
      Automatic mesh
       0
       Monkhorst Pack
       1 1 1
       0. 0. 0.

      POSCAR:
      H2 molecule in a box
      1.0 ! universal scaling parameters
      8.0 0.0 0.0 ! lattice vector a(1)
      0.0 8.0 0.0 ! lattice vector a(2)
      0.0 0.0 8.0 ! lattice vector a(3)
      2 ! number of atoms
      cart ! positions in cartesian coordinates
      0 0 0 ! first atom
      0 0 0.5 ! second atom

      根据定义,E=-(TOTEN(H2)-2*TOTEN(H)) = 6.68eV,键长 = 0.75102 A。比文献中的结合能大了不少。 出现上述结果的原因是在计算单个H原子能量的时候没有指定基态为spin polarized state。为了得到正确的解离能,计算单个H原子能量的时候需要指定ISPIN=2。加上以上的tag后,TOTEN(H atom [spin-polarized])=-1.10351 eV,用公式重新计算解离能:

      E_binding=-(TOTEN(H2)-2*TOTEN(H [spin-polarized ])) = 4.48eV

      和文献吻合。另外,文献[G.Kresse, PRB 62, 8295 (2000)]中详细讨论了氢分子的解离能的计算方法和结果。

    用VASP进行Partial Charge分析实例

    VASP Version : 4.6

    在这篇文章中,我将首先介绍Partial Charge的概念,以及如何用VASP具体的计算Partial Charge。首先,所谓的Partial Charge是针对与Total Charge来说的,指的是某个能量范围、某个K点或者某个特定的态所对应的电荷密度。在文献中最常见的是价带顶部,导带底部,表面态或者局域态所对应的Partial Charge。通过分析这些态所对应的Partial Charge,可以得到体系的一些性质,比如局域态具体的是局域在哪个原子上等。我将通过具体的例子说明如何用VASP进行Partial Charge Analysis。

    进行Partial Charge Analysis的第一步是进行自洽的计算,得到体系的电子结构。这一步的计算采用通常的INCAR和KPOINTS文件。在自洽计算结束后,我们需要保存WAVECAR文件。(通过在INCAR文件中设置LWAVE=TRUE实现)在这个例子中,假设我们需要计算一个硅纳米线的导带和价带的Partial Charge。硅纳米线的结构如下:

     

    http://www.quantumchemistry.net/Experience/UploadFiles_6872/200603/20060331155154521.jpg


    第二步是画出能带结构,以决定你需要画哪条能带的那个K点的态所对应的Partial Charge。关于具体如何用VASP画能带,请参见用VASP4.6计算晶体硅能带实例一文。我们得到硅纳米线的能带结构如下:


    画能带时有些小技巧。你可以用一些支持列模块的编辑器,如UltraEdit,将OUTCAR里的各个K点所对应的本征值粘贴到Origin中。这一步完成后,在Origin中做一个矩阵转置,然后将K点坐标贴到第一列,并将其设为X坐标。如此画出来的基本上就是能带图了。在Origin中可以通过设置纵轴范围来更加清楚的区分费米能级附近的各条能带。如上的硅纳米线所对应的能带结构图如下:

     

    http://www.quantumchemistry.net/Experience/UploadFiles_6872/200603/20060331155540648.jpg
    决定画哪条能带,或者那些感兴趣的K点之后,有如下几种方法计算不同的Partial Charge。如果你希望计算价带顶端的Partial Charge,则需要首先通过能带结构图确定价带的能带标号。需要注意,进行Partial Charge分析必须要保留有自洽计算的WAVECAR才可以。

    第一种Partial Charge分析的INCAR
    ISTART =      1    job   : 0-new  1-cont  2-samecut
    ICHARG =      1   charge: 1-file 2-atom 10-const

    LPARD=.TRUE.
    IBAND= 20 21 22 23
    KPUSE= 1 2 3 4
    LSEPB=.TRUE.
    LSEPK=.TRUE.

    这样的INCAR给出的是指定能带,指定K点所对应的Partial Charge。分析导带、价带等的Partial Charge特性,通常采用的都是这种模式。

    第二种Partial Charge分析的INCAR
    ISTART =      1    job   : 0-new  1-cont  2-samecut
    ICHARG =      1   charge: 1-file 2-atom 10-const

    LPARD=.TRUE.
    EINT = -10.3 -5.1
    LSEPB=.FALSE.
    LSEPK=.FALSE.

    这样的INCAR给出的是在[-10.3 -5.1]能量之间的Partial Charge。这种模式适合于分析某个能量区间内的波函数的性质。

    第三种Partial Charge分析的INCAR
    ISTART =      1    job   : 0-new  1-cont  2-samecut
    ICHARG =      1   charge: 1-file 2-atom 10-const

    LPARD=.TRUE.
    NBMOD=-3
    EINT = -1
    LSEPB=.FALSE.
    LSEPK=.FALSE.

    这样的INCAR给出的是从[Ef-1.0 Ef]能量之间的Partial Charge。这种模式最利于分析费米面附近的波函数的性质。

    用第一种方法,我们可以得到硅纳米线价带顶部和导带底部的Partial Charge如下:

     

    http://www.quantumchemistry.net/Experience/UploadFiles_6872/200603/20060331155753136.jpg