用VASP进行Partial Charge分析实例
用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。
硅纳米线的结构如下：
第二步是画出能带结构，以决定你需要画哪条能带的那个K点的态所对应的Partial Charge。
关于具体如何用VASP画能带，请参见用VASP4.6计算晶体硅能带实例一文。
我们得到硅纳米线的能带结构如下：
画能带时有些小技巧。
你可以用一些支持列模块的编辑器，如UltraEdit，将OUTCAR里的各个K点所对应的本征值粘贴到Origin中。
这一步完成后，在Origin中做一个矩阵转置，然后将K点坐标贴到第一列，并将其设为X坐标。
如此画出来的基本上就是能带图了。
在Origin 中可以通过设置纵轴范围来更加清楚的区分费米能级附近的各条能带。
如上的硅纳米线所对应的能带结构图如下：
决定画哪条能带，或者那些感兴趣的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给出的是在能量之间的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给出的是从能量之间的Partial Charge。
这种模式最利于分析费米面附近的波函数的性质。
用第一种方法，我们可以得到硅纳米线价带顶部和导带底部的Partial Charge如下：
此文的部分结果将发表于Phys. Rev. B。
vasp介绍
vasp介绍(转帖)V ASP (V AMP)主页：http://cms.mpi.univie.ac.at/vasp/简介：V ASP是使用赝势和平面波基组，进行从头量子力学分子动力学计算的软件包，它基于CASTEP 1989版开发。
V AMP/V ASP中的方法基于有限温度下的局域密度近似（用自由能作为变量）以及对每一MD步骤用有效矩阵对角方案和有效Pulay混合求解瞬时电子基态。
这些技术可以避免原始的Car-Parrinello方法存在的一切问题，而后者是基于电子、离子运动方程同时积分的方法。
离子和电子的相互作用超缓V anderbilt赝势(US-PP)或投影扩充波(PA W)方法描述。
两种技术都可以相当程度地减少过渡金属或第一行元素的每个原子所必需的平面波数量。
力与张量可以用V AMP/V ASP很容易地计算，用于把原子衰减到其瞬时基态中。
功能：FeaturesSelf-consistent density functional method with plane wave basisAll-electron projector-augmented-wave (PA W) potentials covering all atoms of the periodic table Local density approximation (LDA) and generalized gradient approximation (GGA)Spin restricted and spin-polarizedSemi-relativistic and full spin-orbit relativisticNon-collinear magnetismLDA (GGA)+U for correlated systemsApplicable to bulk systems, surfaces, interfaces, and molecules (in supercell geometry)Total energies, forces, and the full stress tensorConcurrent relaxation of lattice parameters and atomic positionsAb initio molecular dynamicsGeneration of Monkhorst-Pack special k-pointsIntegration in k-space with smearing or tetrahedron method (with Blöchl correction)Band structure (spin restricted and spin polarized)Site, spin and partial-wave projected densities of statesPartial densities of states at Fermi levelCharge densities Spin densitiesUltra-soft pseudopotentials are offered for comparison平台：UNIX/LINUX相关软件：1. p4vasp主页：http://cms.mpi.univie.ac.at/odubay/p4vasp_site/news.php说明：处理V ASP xml格式输出文件的免费图形环境工具。
VASP计算实例
VASP计算实例目录一、氢气分子H2键长的计算 (3)1.基本文件 (3)2.赝势类型的选择 (3)3.截断能ENCUT参数的选择 (4)4.KPOINTS参数选择 (5)5.对晶格常数进行优化 (6)二、Si晶体晶格常数计算 (8)1.赝势类型选择 (8)2.截断能（ENCUT）参数的选定 (9)3.KPOINTS参数的选定 (11)4.SIGMA参数的选定 (12)5.晶格常数计算结果 (13)三、Si元素单原子能量计算 (14)1.由内聚能倒推单原子能量 (14)2.基本文件 (15)3.单原子能量计算 (15)四、Si的VASP力学常数计算 (16)1.计算所需文件 (16)2.计算与数据处理 (17)3.计算所用到的公式： (18)五、SI晶体的电子结构 (19)1.采用VASP计算能带的步骤 (19)2.电荷分布计算结果 (20)能带计算和结果 (21)3.态密度计算和结果 (21)六、Si晶体介电函数和光学性质的计算 (22)1.计算步骤 (22)2.用到的文件 (23)3.计算结果 (26)七、VASP的声子谱计算 (29)1.计算步骤 (29)2.基本文件 (30)3.声子谱、声子态密度计算和结果 (33)4.热学性质计算和结果 (34)八、化合物co2键长计算 (35)1.计算步骤 (35)2.基本文件 (35)一、氢气分子H2键长的计算1.基本文件准备基本文件INCAR、POTCAR、POSCAR、KPOINT以及脚本文件encut、k、optimize2.赝势类型的选择输入文件如下其中参数要靠经验初选INCAR:System = F2ISTART = 0ICHARG = 2NELMDL = 5ISMEAR = 0SIGMA = 0.1PREC = AccurateKPOINTS：Automatic meshM1 1 10 0 0POSCAR：O115.0 0.00 0.000.00 14.0 0.000.00 0.00 13.01D0.00 0.00 0.00分别选用五个贋势文件进行计算。
vasp经验总结
加U（LMCO为例）
LDAU=.TRUE. LDAUTYPE=2 #The simplified (rotationally invariant) approach to the LSDA+U
#specifies the l-quantum number for which the on-site LDAUL= -1 2 2 -1 interaction is added (-1=no on-site terms added, 1= p, 2= d, 3= f, Default: LDAUL=2)
bnd000.dat文件的第一列数据是k点距离的绝对值，第二列数据是以 Ferim level为参考的本征值。
简立方的高对称点
能量本征值文件EIGENVAL
第一行，前三个整数无意义，第四个整数，如果是2, 表 示是自旋极化的计算，如果是1, 表示非自旋极化的计算。
60 60 1 1 0.1158924E+02 0.5522500E-09 0.2280000E-08 0.5522500E-09 0.5000000E-15 1.000000000000000E-004 第2至5行的数据含义不大明确，可以不管它。 CAR YFTO 492 32 300 第6行的数据表示：第一个数表示体系总的价电子数目，第二个数 表示的计算能带时总的k点数目，第三个数表示的是计算能带时计 算了多少条能带。 0.3742002E-15 0.1262379E-14 0.3742002E-15 0.1388889E-01 1 -46.511682 第8行的前三个数是k点的坐标，第四个数是相 2 -46.511280 应k点的权重。 3 -46.451176 4 -46.449273 第9行给出的是该k点对应的本征值的序号(即第几 5 -46.199455 条能带)，及相应的本征值。 6 -46.197355
VASP用VASP46计算晶体硅能带实例用第一原理计算和其它方法研究纳米体系碳纳米管等的电子和自旋结构
Electronic relaxation 2 (details)
Write flags LWAVE = LCHARG =
T write WAVECAR T write CHGCAR
VASP给INCAR文件中的很多参数都设置了默认值，所以如果你对参 数不熟悉，可以直接用默认的参数值。
vaspviennaabinitiosimulationpackageintroductionoutlinelianghunannormaluniversity用vasp46计算晶体硅能带实例用第一原理计算和其它方法研究纳米体系碳纳米管等的电子和自旋结构参考书记introductionvasp是使用赝势和平面波基组进行从头量子力学分子动力学计算的软件包它基于castep1989版开发
P. 10
KPOINTS 文件：
采用自动的Monkhorst-Pack K点撒取方式。对于类似于硅晶体的
半导体材料，通常 4x4x4 的K点网格就够了。
Monkhorst Pack 0 Monkhorst Pack 44 4 00 0
ISIF = 2
stress and relaxation
P. 11
P. 5
第一行就是K点的倒空间的坐标，接下来的8行告诉我们 在那个K点上的8个能级。你可以通过EXCEL或者ORIGIN 之类的画图软件可视化结果。由于现在手头上已经有了 每个K点的能级信息，则将这些K点的能级连接起来就是 所需要的能带图了。下图是用以上步骤算得的硅的能带 图。我们可以看到硅并非是直接能隙的材料。同时，由 于我们用了LDA，所以硅的能隙和实验相比大大被低估 了（实验为1.12eV，计算值~0.6eV）。
vasp经验总结
POSCAR.sh chmod +x POSCAR.sh 改为可执行文件 ./POSCAR.sh 运行
VASP中，用Berry Phase的方法计算极化值
INCAR中要添加的参数： LCALCPOL=.TRUE. 打开计算Berry的开关 EDIFF=1E-6 精度 DIPOL=0.4 0.4 0.4 选取参考点，任意选的，但是不要和离子重合
BP计算极化过程中，Dipole设置问题
设置在计算离子的dipole时的参考点即设置 DIPOL（注意的是，它的 设置需要使得原子移动前后的原子都在这个参考点的一侧。比如这个例 子中 Al处于(0,0,0)，As处于(0.25, 0.25, 0.25)位置，而将DIPOL设置为( 0.5, 0.5, 0.5)和(0.125, 0.125, 0.125)都是可以的，但是在考虑移动Al原子时，不要将
NELM = 40 # maximum of 40 electronic steps
杂化泛函的计算 --HF Functional
GGA=PS(选用的赝势文件为PBBEsol, =PE为PBE的方法) LHFCALC = .True. PRECFOCK = Normal # NKRED = 2 (设置此参数容易报错，不知为何) TIME = 0.4 HFSCREEN = 0.2 AEXX = 0.25 #the exact exchange is used
2、ALGO, IALGO, LDIAG If the self-consistency loop does not converge within 40 steps, it will probably not converge at all. In this case you should reconsider the tags IALGO, LDIAG, and the mixing-parameters. 一般情况下，或使用IALGO=48时遇到收敛问题的话，可以考虑设IALGO为38， 或设置ALGO=Normal or Fast (in VAS P.4.5 and later versions)。 Default ALGO = Normal 3、NELMDL NELMDL gives the number of non-selfconsistent steps at the beginning In some cases (for instance MD’s, or ionic relaxation) you might set NELMIN to a larger value (4 to 8)
vasp计算
LPARD = .TRUE.开关
IBAND = 480 481 VBM CBM
NBMOD = 1默认
KPUSE = 1第几个K点
LSEPB = .TRUE. vasp查
LSEPK = .TRUE.
20,计算磁性的体系，设置
IALGO=38
LREAL=F
LPLANE=.TRUE.
GGA_COMPAT = F
LMAXMIX=4
LDAU = .TRUE.
LDAUTYPE = 1
LDAUL = 2 -1
LDAUU = 4 0
LDAUJ = 1 0
LORBIT=11
算得更准
19如果要用wannier算磁性体系。Spin要分开，这是应该用vasp+wannier+soc的版本跑，同时INCAR设置成spin=2的模式，但是LSORBIT不开。跑完会出来up和down，再分别对他们局域化。
phonopy-p-c POSCAR-unitcell band.conf
一般来说，对于金属，或者窄能隙半导体，如果用位移法，则需要很大的胞才能算准，但是用DFPT则可以小包算准。对于金属，PBE可能更好点。
9,如果体系较大，EDIFF达到停止计算，很可能是K点取太多，内存不够。
10, bandplot --gnuplot band.yaml >> phon.dat，用origin做声子谱
16,计算真空能级，功函数。在INCAR里面加入LVTOT= T，然后把LOCPOT复制为PARCHG，第一列写成原子名字，第二个原子名字去掉。运行post_VASP，选7即可得到结果。(如果不收敛，可以加大真空层厚度试试)
Partial Charge
用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。
硅纳米线的结构如下：第二步是画出能带结构，以决定你需要画哪条能带的那个K点的态所对应的Partial Charge。
关于具体如何用VASP画能带，请参见用VASP4.6计算晶体硅能带实例一文。
我们得到硅纳米线的能带结构如下：画能带时有些小技巧。
你可以用一些支持列模块的编辑器，如UltraEdit，将OUTCAR里的各个K点所对应的本征值粘贴到Origin中。
这一步完成后，在Origin 中做一个矩阵转置，然后将K点坐标贴到第一列，并将其设为X坐标。
如此画出来的基本上就是能带图了。
在Origin中可以通过设置纵轴范围来更加清楚的区分费米能级附近的各条能带。
如上的硅纳米线所对应的能带结构图如下：决定画哪条能带，或者那些感兴趣的K点之后，有如下几种方法计算不同的Partial Charge。
Vasp入门+实例
0.6 0.5
DOS
0.4 0.3 0.2 0.1 -6 -4 -2 0 2 4 6 8 10
Energy
(4). 做非自洽计算, 求电子结构
• 修改INCAR文件: 将参数ICHARG设为 11 • 修改KPOINTS输入文件
• 运行VASP程序,从输出文件EIGENVAL中提出电子结构
� 1� a1 � a ( i � 2 � 1� a2 � a ( i � 2 � � a3 � ck
第一原理电子结构计算程序：VASP
• 程序原理
• 输入文件
• 输出文件 • 应用
输入文件
POTCAR KPOINTS POSCAR INCAR
Choosing POTCAR file LDA GGA PAW_LDA PAW_GGA PAW_PBE(VASP4.5)
pseudopotentail file Brillouin zone sampling structural data steering parameters
POSCAR输入文件: 原胞中的原子位置
Diamond Si 3.9 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
基矢的公因子
基矢a1 基矢a2
基矢a3 原胞中的原子个数 坐标系选为基矢构成的坐标系
基矢坐标系下原子的位置
� 1 � � a1 � a( j � k ) 2 � 1 � � a2 � a(i � k ) 2 � 1 � � a3 � a(i � j ) 2
1
0
-1
0 0.07 0.14 0.21 0.28 0.34 0.41 0.48 0.55
CHGCAR--vasp手册
CHGCAR £leThis £le contains the lattice vectors, atomic coordinates, the total charge density multiplied by the volume (r) ¤V cell on the £ne FFT-grid (NG(X,Y,Z)F), and the PAW one-center occupancies. CHGCAR can be used to restart VASP from an existing charge density, for visualisation the CHG £le should be used, since the PAW-one centre occupancies are dif£cult to parse. It is possible to avoid that the CHGCAR is written by settingLCHARG = .FALSE.in the INCAR £le (see section 7.45). In VASP, the density is written using the following commands in Fortran:WRITE(IU,FORM) (((C(NX,NY,NZ),NX=1,NGXC),NY=1,NGYZ),NZ=1,NGZC)The x index is the fastest index, and the z index the slowest index. The £le can be read format-free, because at least in new versions, it is guaranteed that spaces separate each number. Please do not forget to divide by the volume before visualizing the £le!For spinpolarized calculations, two sets of data can be found in the CHGCAR £le. The £rst set contains the total charge density (spin up plus spin down), the second one the magnetization density (spin up minus spin down). For non collinear calculations the CHGCAR £le contains the total charge density and the magnetisation density in the x, y and z direction in this order.For dynamic simulation (IBRION=0), the charge density on the £le is the predicted charge density for the next step: i.e. it is compatible with CONTCAR, but incompatible with the last positions in the OUTCAR £le. This allows the CHGCAR and the CONTCAR £le to be used consistently for a molecular dynamics continuation job. For static calculations and relaxations (IBRION=-1,1,2) the written charge density is the selfconsistent charge density for the last step and might be used e.g. for accurate band-structure calculations (see section 10.3).Mind: Since the charge density written to the £le CHGCAR is not the selfconsistent chargedensity for the positions on the CONTCAR £le, do not perform a bandstructure calculation (ICHARG=11) directly after a dynamic simulation (IBRION=0) (see section 10.3).Band decomposed charge density (parameters) 分波电荷密度计算VASP can calculate the partial (band decomposed) charge density部分（分波）电荷密度 according to parameters specified in the file INCAR. It must be noted, that the densities calculated by VASP (including the band decomposed charge density) are always symmetrized对称的 using the space group空间群 and point group symmetry determined by VASP.根据设置INCAR中的参数，VASP可以计算分波电荷密度，算分波电荷密度通常考虑空间群和点群对称性。
Partial charge density分析(实操版)
一、含义Partial charge density：也称为Band decomposed charge densityPartial charge density计算，就是计算特定的某个(或某些)k点和本征值(这些k点和本征值是相互对应的)所对应的本征波函数的平方(也就是电荷密度)。
二、计算过程和参数设置计算Partial电荷密度的方法：1.结构驰豫2.计算静态，输出CHGCAR、W A VECAR3.计算能带，无须输出CHGCAR、W A VECAR4.做能带图，选点。
5.计算Partial电荷密度，在INCAR中设置需要计算的能带编号和K点编号，不需要KPOINTS文件（也可以放一个算静态或者算能带的KPOINTS文件，都不影响结果），读入静态的CHGCAR和WA VECAR。
6.计算结束后，用PARCHG.****.****文件删掉后面数字，改成PARCHG作图。
参数设置：以上第一步和以往的计算一样，参数就不介绍了。
第4步：选点1.需要一个小程序（band_plot来处理能带计算输出的OUTCAR和EIGENV AL），或者一个可执行脚本（band_plot放到计算目录的bin/目录下），两者可以相互转化。
此处，用后面这种，将band_plot放到./bin目录下：2.然后到能带计算完成的目录下，输入命令：band_plot3.然后会得到一个文件数据文件：band_plot.dat4.将这个文件下载，用origin作图5.调整图片坐标6.我们根据需要在图中选择一点，然后找出其坐标，为了操作方便需要放大图片View—Zoom In:为了选点方便，采用点模式用左侧快捷键来得到需要点的坐标，如图【X=04714，Y=0.41041】然后，我们到数据表格中，搜索0.41041这个数，（因为X=04714是第一列，一目了然）Edit-Find，出现对话框，输入需要找的数：从表中，我们得知，第一列，对应X轴，也对应K点，我们算了31个K点，我们找的那一点，正好是第31个K点。