http://abhiramn.blogspot.com/2010/06/latex-on-blogspot.html
$\Gamma$
$ \mathcal{E}_\mathrm{defect} = E_\mathrm{defect}- E_\mathrm{host}+n_\mathrm{C}\mu_\mathrm{C}-n_\mathrm{H}\mu_\mathrm{H}, $
Friday, September 30, 2011
Wednesday, June 18, 2008
如何用vasp计算Born有效电荷
VASP是采用的Berry phase方法来计算半导体和绝缘体材料的Born有效电荷。这种方法的介绍可以参考VASP手册上提到的文献。
如VASP手册上的介绍，在采用Berry phase方法来计算某个原子某个方向的Born有效电荷时，有两大步：
1)计算离子未移动时(即平衡状态时)的极化(polarization); 2)手动移动所要计算的原子的位置，计算由原子偏离平衡态时的极化。每一大步中有包括了四个小步骤: i)自恰计算体系以得到收敛的电荷密度(CHG, CHGCAR文件)和波函数(WAVECAR); ii)计算G1方向上的电子极化; iii)计算G2方向上的电子极化; iv)计算G3方向上的电子极化。
最后从每一小步中的OUTCAR文件中找出极化各个部分的贡献:ev (电子项), bp (Berry phase项), ion (离子项)，进行简单的数值加减计算，按差分方法计算得到Born有效电荷。 对每类原子的各个方向按同样的步骤进行计算处理，在计算时可根据晶体的对称性，只需要计算Born有效电荷张量的某些分量即可。
下面以计算AlAs的As原子的Born有效电荷为例来介绍具体操作步骤。
1、计算AlAs平衡态时的极化：
i)自恰计算平衡态时的AlAs，主要的输入文件：
INCAR:
SYSTEM = AlAs
ENCUT = 300
ISTART = 0
ICHARG = 2
ISMEAR = 0; SIGMA = 0.002
EDIFF = 1E-5
PREC = Accurate
(注释: 这里采用ISMEAR=0; SIGMA=0.002，是因为所计算的体系是半导体或绝缘体，电子的占有数实际上固定的，非0就是2或1了，通过将ISMAR=0(即Gaussian方法来确定电子占有数)且SIGMA设置为一个很小的值来达到一种类似PWSCF程序中的occupation-fixed 方法的效果。如果参看VASP的源代码，可以看到：代码里面会根据INCAR中设置了LBERRY = .TRUE.而自动将ISMEAR和SIGMA进行调整为默认的值:ISMEAR = 0; SIGMA = 0.00)。
KPOINTS:
4x4x4
0
Monma
4 4 4
0 0 0
POSCAR:
AlAs
5.57223
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1 1
Direct
0.000 0.0000 0.0000
0.250 0.2500 0.2500
POTCAR，采用的是LDA的ultrasoft pseudopotential。
自恰计算完之后将得到的CHGCAR, CHG和WAVECAR保存。
ii)计算沿着G1方向上的电子极化：
输入文件如下：
INCAR:
SYSTEM = AlAs
ENCUT = 300
#ISTART = 0
#ICHARG = 2
ISMEAR = 0; SIGMA = 0.002
EDIFF = 1E-5
PREC = Accurate
LBERRY = .TRUE.
IGPAR = 1
NPPSTR = 8
DIPOL = 0.5 0.5 0.5
（注释：要采用Berry phase方法来计算电子极化时，需在INCAR中设置LBERRY=.TRUE.，另外设置所要计算的方向即设置IGPAR （可赋的值为1,2,3分别表示G1, G2和G3方向), 沿着IGPAR方向上的一串k点的个数 即设置NPPSTR(注意的是，在平衡态时的计算和原子移动之后的计算中，它值应该是一样的)，最后设置在计算离子的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原子时，不要将原子移动原胞之外即偏移量为负数；另外也不要将DIPOL设置在所要移动的原子上，如果是这样的话，则会导致移动该原子后，该原子不在DIPOL的同一侧，使得原子移动之后的极化的Berry-phase项要比平衡态时的大很多。）
POSCAR和KPOINTS同i)中的。
将第i)计算得到的CHGCAR, CHG, WAVECAR拷贝过来，进行计算。
计算完后保存OUTCAR，从OUTCAR文件中找出极化的各个项的贡献，可以采用下面的命令：
grep '' OUTCAR
如下面的数据
Expectation value term:《R》ev
《R》 = ( -0.00022, -0.00117, 0.00074 ) electrons Angst
Berry-Phase term: bp
= ( 0.00000, -0.00455, -0.00455 ) electrons Angst
ionic term: ion
= ( 15.32363, 15.32363, 15.32363 ) electrons Angst
(注释：后面的三列数据分别表示的是x, y, z方向的。)
iii)和iv)步同ii)类似，只是分别将IGPAR = 1改为IGPAR = 2和IGPAR = 3。
下面是计算G2方向时(即IGPAR=2)得到的数据:
Expectation value term: ev
= ( -0.00117, -0.00022, 0.00074 ) electrons Angst
Berry-Phase term: bp
= ( -0.00455, 0.00000, -0.00455 ) electrons Angst
ionic term: ion
= ( 15.32363, 15.32363, 15.32363 ) electrons Angst
下面是计算G3方向时(即IGPAR=3)得到的数据：
Expectation value term: ev
= ( -0.00117, 0.00074, -0.00022 ) electrons Angst
Berry-Phase term: bp
= ( -0.00455, -0.00455, 0.00000 ) electrons Angst
ionic term: ion
= ( 15.32363, 15.32363, 15.32363 ) electrons Angst
下面按公式计算计平衡态时 ev在G1、G2和G3方向上的平均值：
ev, undistorted, average =
[( -0.00022, -0.00117, 0.00074 )+
( -0.00117, -0.00022, 0.00074 )+
( -0.00117, 0.00074, -0.00022 ) ]/3.0
= ( -0.00022-0.00117-0.00117, -0.00117-0.00022+0.00074, 0.00074+0.00074-0.00022 ) /3.0
=(-0.000853333, -0.000216667, 0.00042)
把平衡态时各个方向计算时得到 bp项求和：
bp, undistorted, all =
( 0.00000, -0.00455, -0.00455 ) +
( -0.00455, 0.00000, -0.00455 ) +
( -0.00455, -0.00455, 0.00000 )
= ( -0.00455*2, -0.00455*2, -0.00455*2 )
= ( -0.0091, -0.0091, -0.0091)
那么平衡态时电子极化项 el, undistorted就是:
el, undistorted = ev, undistorted, average + bp, undistorted, all
=(-0.000853333, -0.000216667, 0.00042) +( -0.0091, -0.0091, -0.0091)
=(-0.000853333 -0.0091, -0.000216667 -0.0091, 0.00042 -0.0091)
=(-0.00995333, -0.00931667, -0.00868)
平衡态时离子极化项 ion, undistorted，在任意一个方向计算时得到的数据，即：
ion, undistorted = ( 15.32363, 15.32363, 15.32363 )
2. 计算Al原子偏离平衡位置时的极化：
i)其中的INCAR， KPOINTS同第1 i)步中的，只是将POSCAR中 As原子的位置沿着某个方向略移动一个小量，如下：
POSCAR
5.57223
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1 1
Direct
0.000 0.0000 0.0000
0.250 0.2500 0.2300
（注释：这种设置是导致As在x和y方向上都移动了0.05572 angstrom，这个可以自己算出来，或者直接从移动前后两次计算的OUTCAR中的原子cartesian坐标进行简单的加减计算之后得到。)
然后作自恰计算，计算之后保存得到CHGCAR, CHG, WAVECAR。
ii)、计算原子移动之后，体系沿着G1方向的极化：
INCAR和KPOINTS的输入文件同第1 ii)步中的，将上一步保存下来的CHG, CHCAR和WAVECAR拷贝过来，然后计算得到OUTCAR，在保存后，从中可以找到相关的极化计算数据，通过下面的命令来查找：
grep '' OUTCAR
此例子的计算数据：
Expectation value term: 《R》ev
《R》 = ( -0.00091, -0.00100, -0.00030 ) electrons Angst
Berry-Phase term: 《R》bp
《R》 = ( 0.00000, -0.00475, -0.00475 ) electrons Angst
ionic term: ion
《R》= ( 15.60224, 15.60224, 15.32363 ) electrons Angst
iii)& iv).计算原子移动后体系沿着G2和G3方向的极化：
INCAR和KPOINTS中的分别与第1 iii)&iv)步中的一样。
此例子计算得到的数据:
沿着G2方向的：
Expectation value term: 《R》ev
《R》 = ( -0.00100, -0.00091, -0.00030 ) electrons Angst
Berry-Phase term: bp
《R》= ( -0.00475, 0.00000, -0.00475 ) electrons Angst
ionic term: ion
《R》 = ( 15.60224, 15.60224, 15.32363 ) electrons Angst
沿着G3方向的：
Expectation value term: 《R》ev
《R》 = ( -0.00233, 0.00032, -0.00029 ) electrons Angst
Berry-Phase term: 《R》bp
《R》 = ( -0.40026, -0.40026, 0.00000 ) electrons Angst
ionic term: 《R》ion
《R》 = ( 15.60224, 15.60224, 15.32363 ) electrons Angst
同样按公式计算出原子移动之后的极化：
将G1, G2, G3各个方向计算得到 《R》ev项进行平均：
《R》ev, distorted, average =
[ ( -0.00091, -0.00100, -0.00030 ) +
( -0.00100, -0.00091, -0.00030 ) +
( -0.00233, 0.00032, -0.00029 )]/3.0
= (
=(-0.00141333, -0.00053, -0.000103333)
将G1, G2, G3各个方向计算时得到《R》bp项求和：
《R》bp, distorted, all =
( 0.00000, -0.00475, -0.00475 )+
( -0.00475, 0.00000, -0.00475 )+
( -0.40026, -0.40026, 0.00000 )
=(0.00000-0.00475-0.40026, -0.00475+0.00000-0.40026, -0.00475 -0.00475+0.00000)
=(-0.40501, -0.40501, -0.0095)
由此得到原子移动之后电子极化项 《R》el, distorted:
《R》el, distorted= 《R》ev, distorted, average + 《R》bp, distorted, all
=(-0.00141333, -0.00053, -0.000103333)+(-0.40501, -0.40501, -0.0095)
=(-0.00141333-0.40501, -0.00053-0.40501, -0.000103333-0.0095)
=(-0.406423, -0.40554, -0.00960333)
原子移动后离子极化项《R》ion, distorted取任意一个方向上计算得到的数据《R》ion:
《R》ion, distorted = ( 15.60224, 15.60224, 15.32363 )
由此得到由于As原子的移动引起的极化变化为：
Delta 《R》= 《R》ion, distorted + 《R》el, distorted - 《R》ion, undistorted - 《R》el, undistorted
= ( 15.60224, 15.60224, 15.32363 )+(-0.406423, -0.40554, -0.00960333) - ( 15.32363, 15.32363, 15.32363 )-(-0.00995333, -0.00931667, -0.00868)
=( 15.60224-0.406423-15.32363+0.00995333, 15.60224-0.40554-15.32363+0.00931667, 15.32363-0.00960333-15.32363+0.00868)
=(-0.11786, -0.117613, -0.00092333)
在前面已经提到了这个是在As分别在沿着x, y方向上移动了0.05572 angstrom时的计算结果。那么由此可计算出As原子的沿着x,y方向上的Born有效电荷（实际上由于AlAs是立方晶体，As原子的Born有效电荷沿着xx, yy, zz方向的都相等，xz, yz, xy方向上的都为0）按差分公式可以计算得到：
Z*_As = Delta 《R》 / Delta x_As
=-0.11786/0.05572
=-2.11522
由于原胞所有原子的Born有效电荷各个方向的加起来要满足电荷中性的条件，AlAs中只有Al和As两类(也只有2个原子)原子，因此也可以得到Al的Born有效电荷为2.11522，这与ABINIT中的教程中计算的数据很接近。
实际上面两大步的4个小步骤，可采用下面的脚本文件来连续计算完成：
要注意的就是DIPOL 的设置，以及原子移动的选取：不要将DIPOL设置在某个原子位置上，原子移动的选取要保证原子在移动前后都是DIPOL的同一侧。
如VASP手册上的介绍，在采用Berry phase方法来计算某个原子某个方向的Born有效电荷时，有两大步：
1)计算离子未移动时(即平衡状态时)的极化(polarization); 2)手动移动所要计算的原子的位置，计算由原子偏离平衡态时的极化。每一大步中有包括了四个小步骤: i)自恰计算体系以得到收敛的电荷密度(CHG, CHGCAR文件)和波函数(WAVECAR); ii)计算G1方向上的电子极化; iii)计算G2方向上的电子极化; iv)计算G3方向上的电子极化。
最后从每一小步中的OUTCAR文件中找出极化各个部分的贡献:
下面以计算AlAs的As原子的Born有效电荷为例来介绍具体操作步骤。
1、计算AlAs平衡态时的极化：
i)自恰计算平衡态时的AlAs，主要的输入文件：
INCAR:
SYSTEM = AlAs
ENCUT = 300
ISTART = 0
ICHARG = 2
ISMEAR = 0; SIGMA = 0.002
EDIFF = 1E-5
PREC = Accurate
(注释: 这里采用ISMEAR=0; SIGMA=0.002，是因为所计算的体系是半导体或绝缘体，电子的占有数实际上固定的，非0就是2或1了，通过将ISMAR=0(即Gaussian方法来确定电子占有数)且SIGMA设置为一个很小的值来达到一种类似PWSCF程序中的occupation-fixed 方法的效果。如果参看VASP的源代码，可以看到：代码里面会根据INCAR中设置了LBERRY = .TRUE.而自动将ISMEAR和SIGMA进行调整为默认的值:ISMEAR = 0; SIGMA = 0.00)。
KPOINTS:
4x4x4
0
Monma
4 4 4
0 0 0
POSCAR:
AlAs
5.57223
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1 1
Direct
0.000 0.0000 0.0000
0.250 0.2500 0.2500
POTCAR，采用的是LDA的ultrasoft pseudopotential。
自恰计算完之后将得到的CHGCAR, CHG和WAVECAR保存。
ii)计算沿着G1方向上的电子极化：
输入文件如下：
INCAR:
SYSTEM = AlAs
ENCUT = 300
#ISTART = 0
#ICHARG = 2
ISMEAR = 0; SIGMA = 0.002
EDIFF = 1E-5
PREC = Accurate
LBERRY = .TRUE.
IGPAR = 1
NPPSTR = 8
DIPOL = 0.5 0.5 0.5
（注释：要采用Berry phase方法来计算电子极化时，需在INCAR中设置LBERRY=.TRUE.，另外设置所要计算的方向即设置IGPAR （可赋的值为1,2,3分别表示G1, G2和G3方向), 沿着IGPAR方向上的一串k点的个数 即设置NPPSTR(注意的是，在平衡态时的计算和原子移动之后的计算中，它值应该是一样的)，最后设置在计算离子的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原子时，不要将原子移动原胞之外即偏移量为负数；另外也不要将DIPOL设置在所要移动的原子上，如果是这样的话，则会导致移动该原子后，该原子不在DIPOL的同一侧，使得原子移动之后的极化的Berry-phase项要比平衡态时的大很多。）
POSCAR和KPOINTS同i)中的。
将第i)计算得到的CHGCAR, CHG, WAVECAR拷贝过来，进行计算。
计算完后保存OUTCAR，从OUTCAR文件中找出极化
grep '
如下面的数据
Expectation value term:
Berry-Phase term:
ionic term:
(注释：
iii)和iv)步同ii)类似，只是分别将IGPAR = 1改为IGPAR = 2和IGPAR = 3。
下面是计算G2方向时(即IGPAR=2)得到的数据:
Expectation value term:
Berry-Phase term:
ionic term:
下面是计算G3方向时(即IGPAR=3)得到的数据：
Expectation value term:
Berry-Phase term:
ionic term:
下面按公式计算计平衡态时
[( -0.00022, -0.00117, 0.00074 )+
( -0.00117, -0.00022, 0.00074 )+
( -0.00117, 0.00074, -0.00022 ) ]/3.0
= ( -0.00022-0.00117-0.00117, -0.00117-0.00022+0.00074, 0.00074+0.00074-0.00022 ) /3.0
=(-0.000853333, -0.000216667, 0.00042)
把平衡态时各个方向计算时得到
( 0.00000, -0.00455, -0.00455 ) +
( -0.00455, 0.00000, -0.00455 ) +
( -0.00455, -0.00455, 0.00000 )
= ( -0.00455*2, -0.00455*2, -0.00455*2 )
= ( -0.0091, -0.0091, -0.0091)
那么平衡态时电子极化项
=(-0.000853333, -0.000216667, 0.00042) +( -0.0091, -0.0091, -0.0091)
=(-0.000853333 -0.0091, -0.000216667 -0.0091, 0.00042 -0.0091)
=(-0.00995333, -0.00931667, -0.00868)
平衡态时离子极化项
2. 计算Al原子偏离平衡位置时的极化：
i)其中的INCAR， KPOINTS同第1 i)步中的，只是将POSCAR中 As原子的位置沿着某个方向略移动一个小量，如下：
POSCAR
5.57223
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1 1
Direct
0.000 0.0000 0.0000
0.250 0.2500 0.2300
（注释：这种设置是导致As在x和y方向上都移动了0.05572 angstrom，这个可以自己算出来，或者直接从移动前后两次计算的OUTCAR中的原子cartesian坐标进行简单的加减计算之后得到。)
然后作自恰计算，计算之后保存得到CHGCAR, CHG, WAVECAR。
ii)、计算原子移动之后，体系沿着G1方向的极化：
INCAR和KPOINTS的输入文件同第1 ii)步中的，将上一步保存下来的CHG, CHCAR和WAVECAR拷贝过来，然后计算得到OUTCAR，在保存后，从中可以找到相关的极化计算数据，通过下面的命令来查找：
grep '
此例子的计算数据：
Expectation value term:
Berry-Phase term:
ionic term:
iii)& iv).计算原子移动后体系沿着G2和G3方向的极化：
INCAR和KPOINTS中的分别与第1 iii)&iv)步中的一样。
此例子计算得到的数据:
沿着G2方向的：
Expectation value term:
Berry-Phase term:
ionic term:
沿着G3方向的：
Expectation value term:
Berry-Phase term:
ionic term:
同样按公式计算出原子移动之后的极化：
将G1, G2, G3各个方向计算得到
[
( -0.00100, -0.00091, -0.00030 ) +
( -0.00233, 0.00032, -0.00029 )]/3.0
= (
=(-0.00141333, -0.00053, -0.000103333)
将G1, G2, G3各个方向计算时得到《R》bp项求和：
《R》bp, distorted, all =
( 0.00000, -0.00475, -0.00475 )+
( -0.00475, 0.00000, -0.00475 )+
( -0.40026, -0.40026, 0.00000 )
=(0.00000-0.00475-0.40026, -0.00475+0.00000-0.40026, -0.00475 -0.00475+0.00000)
=(-0.40501, -0.40501, -0.0095)
由此得到原子移动之后电子极化项 《R》el, distorted:
《R》el, distorted= 《R》ev, distorted, average + 《R》bp, distorted, all
=(-0.00141333, -0.00053, -0.000103333)+(-0.40501, -0.40501, -0.0095)
=(-0.00141333-0.40501, -0.00053-0.40501, -0.000103333-0.0095)
=(-0.406423, -0.40554, -0.00960333)
原子移动后离子极化项《R》ion, distorted取任意一个方向上计算得到的数据《R》ion:
《R》ion, distorted = ( 15.60224, 15.60224, 15.32363 )
由此得到由于As原子的移动引起的极化变化为：
Delta 《R》= 《R》ion, distorted + 《R》el, distorted - 《R》ion, undistorted - 《R》el, undistorted
= ( 15.60224, 15.60224, 15.32363 )+(-0.406423, -0.40554, -0.00960333) - ( 15.32363, 15.32363, 15.32363 )-(-0.00995333, -0.00931667, -0.00868)
=( 15.60224-0.406423-15.32363+0.00995333, 15.60224-0.40554-15.32363+0.00931667, 15.32363-0.00960333-15.32363+0.00868)
=(-0.11786, -0.117613, -0.00092333)
在前面已经提到了这个是在As分别在沿着x, y方向上移动了0.05572 angstrom时的计算结果。那么由此可计算出As原子的沿着x,y方向上的Born有效电荷（实际上由于AlAs是立方晶体，As原子的Born有效电荷沿着xx, yy, zz方向的都相等，xz, yz, xy方向上的都为0）按差分公式可以计算得到：
Z*_As = Delta 《R》 / Delta x_As
=-0.11786/0.05572
=-2.11522
由于原胞所有原子的Born有效电荷各个方向的加起来要满足电荷中性的条件，AlAs中只有Al和As两类(也只有2个原子)原子，因此也可以得到Al的Born有效电荷为2.11522，这与ABINIT中的教程中计算的数据很接近。
实际上面两大步的4个小步骤，可采用下面的脚本文件来连续计算完成：
要注意的就是DIPOL 的设置，以及原子移动的选取：不要将DIPOL设置在某个原子位置上，原子移动的选取要保证原子在移动前后都是DIPOL的同一侧。
Subscribe to:
Posts (Atom)