4.2 光学性质

介电常数

首先进行结构优化,参考官网: https://www.vasp.at/wiki/index.php/Dielectric_properties_of_SiC

静态介电常数:

方法1:DFPT方法;

{.line-numbers}
1
2
3
4
5
IBRION = 8
NSW=1
LEPSILON = .TRUE.
LRPA=.TRUE.#RPA考虑局域场效应,默认关闭,即IP近似
LPEAD=.TRUE.#有限差分法

方法2:Response to finite electric fields

{.line-numbers}
1
2
3
4
5
LCALCEPS = .TRUE.
IBRION=6
NFREE=2
ISMEAR = 0
SIGMA = 0.01

在OUTCAR中抓取:

STATIC DIELECTRIC TENSOR (including local field effects in DFT)``` 电子贡献或者
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
```MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))```
```MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC CONTRIBUTION``` 离子贡献
```BORN EFFECTIVE CHARGES (in e, cummulative output)```
```ELASTIC MODULI IONIC CONTR (kBar)```
```PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)```
### 频率依赖介电函数:读取WAVECAR
IPA:
```javascript {.line-numbers}
ALGO = Exact
NBANDS = 64
LOPTICS = .TRUE.
CSHIFT = 0.100
NEDOS = 2000
#and you might try with the following
#LPEAD = .TRUE.

frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects)
and
frequency dependent REAL DIELECTRIC FUNCTION (independent particle, no local field effects)

Note:对于光学性质的计算,也就是计算材料的介电函数,需要足够多的空带和致密的K网格点,使其达到非常好的收敛状态,我们才可以得到合理的光学性质;因此通常计算中,一般设置NBANDS为默认值(可在自洽的OUTCAR中以关键字NBANDS查找到)的2-3倍,K网格点为自洽值或适当增加。
使用optics.sh脚本得到介电函数的实部和虚部

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
# extract image and real parts of dielectric function from vasprun.xml
awk 'BEGIN{i=1} /<imag>/,\
/<\/imag>/ \
{a[i]=$2 ; b[i]=$3 ; c[i]=$4; d[i]=$5 ; e[i]=$6 ; f[i]=$7; g[i]=$8; i=i+1} \
END{for (j=12;j<i-3;j++) print a[j],b[j],c[j],d[j],e[j],f[j],g[j]}' vasprun.xml > IMAG.in
#awk 'BEGIN{i=1} /imag/,\
# /\/imag/ \
# {a[i]=$2 ; b[i]=$3 ; c[i]=$4; d[i]=$5 ; e[i]=$6 ; f[i]=$7; g[i]=$8; i=i+1} \
# END{for (j=12;j<i-3;j++) print a[j],b[j],c[j],d[j],e[j],f[j],g[j]}' vasprun.xml > IMAG.in
awk 'BEGIN{i=1} /<real>/,\
/<\/real>/ \
{a[i]=$2 ; b[i]=$3 ; c[i]=$4; d[i]=$5 ; e[i]=$6 ; f[i]=$7; g[i]=$8; i=i+1} \
END{for (j=12;j<i-3;j++) print a[j],b[j],c[j],d[j],e[j],f[j],g[j]}' vasprun.xml > REAL.in

vaspkit处理:
vaspkit-task 711

压电常数

上面计算完介电常数后,一般正交晶系。
grep "PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)" OUTCAR -A10
&
grep "PIEZOELECTRIC TENSOR for field in x, y, z (C/m^2)" OUTCAR -A10

极化强度

  1. 极化强度计算
    首先选取极化相FE(PhaseB, P, 铁电相)和参考相NP(PhaseA, P=0, 中心对称相),分别以下列INCAR计算,
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
     SYSTEM = FE
    ENCUT = 600 #铁电极化计算建议设置高精度
    EDIFF = 1E-6
    IBRION = -1
    POTIM = 0.3
    NSW = 0
    EDIFFG = -5E-3
    ISMEAR = 0
    SIGMA = 0.05
    PREC = Accurate
    ISIF = 2
    #NPAR = 4
    LWAVE = .FALSE.
    LCHARG = .FALSE.
    LREAL = .FALSE.
    LDIPOL = .T.
    DIPOL = 0.2 0.2 0.2 #不要设置在原子及迁移路径上,设置在真空层的一边/质心
    LCALCPOL = .TRUE. #计算铁电极化开关
    IPEAD=4
    LPEAD=.True.
  2. 收集结果:grep “Total electronic dipole momen” OUTCAR & grep “Ionic dipole moment” OUTCAR,离子与电子相加即可,然后用NP相-铁电相即可得到极化强度P,注意单位换算。
    首先要理清数据单位。VASP 计算得到的 dipole moment 的单位是 e*Å,它与库仑之间换算为:
    $1e = 1.602176634 * 10^{-19} C$
    $1Å = 10^{-10} m$
    三维体系的极化强度: 极化值除以体积。单位为 $e/Å^2$;二维体系的极化强度:极化值除以面积。单位为 e/Å。

Note:

在铁电极化计算过程中经常会出现参考相为非半导体的情况,这种情况下可以:①镜像法:以FE以NP为参照中心,建立一个-FE相(PhaseB’, 即铁电极化方向相反)。然后按照上述方法计算极化,再以FE相的极化-(-FE相)的极化值,然后除以2即得到极化强度。可以理解为1-(-1)=2,2/2=1;②线性插值法:在FE和NP相中间插入一系列中间过渡相(0%(FE),10%,20%,30%….100%(NP)),计算它们的极化,然后可以用Origin做拟合,100%时即为NP相的极化。(100%-0%)则是极化强度P。

线性插值法计算极化强度

原文链接:

https://www.cnblogs.com/ghzhan/articles/16305679.html

BaTiO3 是钙钛矿结构,它的铁电相结构和中心对称结构如图所示0401,属于四方晶系。

  1. 先构造 BaTiO3 两种极化方向的晶格结构,并用 VASP 进行结构优化得到 CONTCAR;

  2. 将上一步优化后的两个结构分别放入创建好的 ini, fin 文件夹。利用 NEB 的 nebmake.pl 命令产生这两种极化方向的中间过渡结构 (vtst下载地址: Download — Transition State Tools for VASP (utexas.edu)),具体命令为:

    nebmake.pl ini/CONTCAR fin/CONTCAR 32

    这里的 ‘32’ 是表示产生中间过渡的 32 种结构。执行上述命令后,当前文件夹下会产生 00, 01, 02, …, 33 个文件夹,每个文件夹下有一个 POSCAR 文件。

  3. 对每个文件夹的结构进行一次 VASP 自洽运算,INCAR 文件里面需要额外设置 DIPOL 和 LCALPOL 参数。

    DIPOL 参数可以选取任一坐标,但是保证同一体系采用相同值。

    LCALPOL 参数是打开 极化运算。

    INCAR 文件:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    SYSTEM=BaTiO3
    ISTART =0
    ICHARG =2
    PREC =Accurate
    ENCUT =520
    EDIFF =0.1E-07
    ISMEAR = 0
    LWAVE=.FALSE.
    LCHARG = .FALSE.
    NELM = 200
    DIPOL = 0.5 0.5 0.5
    LCALCPOL=.TRUE
  4. 计算完成后运行命令:

    1
    2
    3
    grep 'dipole moment' OUTCAR|rev| awk '{printf ("%s ", $4)}'|rev|awk '{printf ("%f\n", $1+$2)}' >> ../dipole_c.dat

    grep 'free energy TOTEN' OUTCAR | awk '{printf ("%s\n",$5)}' >> ../energy.dat

    这里通过 grep 命令产生的 dipole_c.dat 文件记录的是沿着 c 方向的极化值,这是因为 BaTiO3 是沿着 c 方向极化的。对于具体的情况需要自行修改。

  5. 处理和分析数据。

    首先要理清数据单位。VASP 计算得到的 dipole moment 的单位是 e*Å,它与库仑之间换算为:

    $1e = 1.602176634 * 10^{-19} C$

    $1Å = 10^{-10} m$

    三维体系的极化强度: 极化值除以体积。单位为 e/Å2;二维体系的极化强度:极化值除以面积。单位为 e/Å。

    在这个例子中,BaTiO3 的体积为 64.3586 Å3。接下来,我们来处理数据。用energy.dat 文件 (单位为 eV),这样,我们就得到每个中间结构 (image) 的总能:0402
    用dipole_c.dat 文件 (单位: e*Å):不同 image 的极化值如图所示0403
    可以发现,该极化强度不是连续的,这是和选取的原胞有关,需要考虑极化量子的影响。BaTiO3 沿着 c 方向极化,所以需要对该极化值加减整数倍的 c 方向的晶格常数。这样我们得到下图:0404

    选取最靠近中间极化强度为 0 的那条曲线,即为 BaTiO3 的极化强度曲线:0405

    除以体积并进行简单的单位换算后为:0406
    另外,我们也可以绘制极化值 P 与能量 E 曲线0407
    因此, BaTiO3 的极化强度大约为 0.248 C/m2,与文献中的实验结果 0.26 C/m2 吻合。(Physical Review, 99(4), 1161–1165, 1955 )

    拟合 Landau-Ginzburg 公式

    $E= \sum_{i} A/2 * P_i^2 + B/4 * P_i^4 + C/6 * P_i^6 + D/2 * \sum_{i,j}(P_i-P_j)^2$

拟合该多项式曲线,得到

$E = 6062.434883706029P^6 + 2437.756598493882P^4 + -340.5806845162749*P^2 + 10.339600058315137$
其中参数 A, B, C 分别为

$A = -0.68116 eV *(m^2/C)^2
B = 9.75103 eV *(m^2/C)^4
C = 36.37461 eV *(m^2/C)^6$
0408

  1. 绘制能量曲线的脚本
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    with open('energy.dat','r') as f:
    content = f.readlines()
    res = [float(i.strip('\n')) for i in content]
    minres = min(res)
    res = [1e3*(i-minres) for i in res]
    import matplotlib.pyplot as plt
    plt.figure()
    plt.plot(res,'b.-')
    plt.xlabel('displacement')
    plt.ylabel('Energy ($meV$)')
    plt.show()

绘制dipole曲线,以及考虑极化量子的脚本:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
with open('dipole_c.dat','r') as f:
content = f.readlines()
res = [float(i.strip('\n')) for i in content]

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('displacement')
plt.ylabel('Polarization ($e*{\AA}$)')
plt.xlim([0,32])
plt.plot(res,'r.')


### 考虑极化量子
c = 4.0335000000000001
tmp = []
N = 20
for j in range(N):
start = -N/2+j
tmp1 = [i+start*c for i in res]
# print(tmp1[0])
tmp.append(tmp1)
plt.figure()
for i in tmp:
plt.plot(i,'.')

#plt.ylim([-2,2])
plt.xlabel('displacement')
plt.ylabel('Polarization ($e*{\AA}$)')
plt.xlim([0,32])

### 单位换算
P = [-1.1965400000000006, -1.1288099999999996, -1.060039999999999, -0.990219999999999, -0.9193599999999993, -0.8474299999999992, -0.7744599999999995, -0.7004900000000003, -0.6255299999999995, -0.5496499999999997, -0.4728999999999992, -0.3953799999999994, -0.31716999999999906, -0.2384000000000004, -0.1591799999999992,-0.11827, 0.0, 0.11827, 0.1591799999999992, 0.2384000000000004, 0.31716999999999906, 0.3953799999999994, 0.4728999999999992, 0.5496499999999997, 0.6255299999999995, 0.7004900000000003, 0.7744599999999995, 0.8474299999999992, 0.9193599999999993, 0.990219999999999, 1.060039999999999, 1.1288099999999996, 1.1965400000000006]
plt.figure()
P = [i*0.248945227832799346 for i in P]
plt.xlabel('displacement')
plt.ylabel('Polarization ($C/m^2$)')
plt.xlim([0,32])
plt.plot(P,'r.')
plt.show()

绘制极化值 P 与能量 E 曲线的脚本:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
E = [0.0035557599999975764, 0.0016252999999935014, 0.0004821100000000911, 0.0, 5.983999999870093e-05, 0.0005505999999968481, 0.0013690299999993272, 0.0024182999999950994, 0.003608839999998281, 0.004859239999994713, 0.006096989999996083, 0.007258279999994954, 0.008287759999994648, 0.009138139999997463, 0.00977216999999797, 0.010162939999993625, 0.010294929999993485, 0.010162939999993625, 0.00977216999999797, 0.009138139999997463, 0.008287759999994648, 0.007258279999994954, 0.006096989999996083, 0.004859239999994713, 0.003608839999998281, 0.0024182999999950994, 0.0013690299999993272, 0.0005505999999968481, 5.983999999870093e-05, 0.0, 0.0004821100000000911, 0.0016252999999935014, 0.0035557599999975764]
P = [-1.1965400000000006, -1.1288099999999996, -1.060039999999999, -0.990219999999999, -0.9193599999999993, -0.8474299999999992, -0.7744599999999995, -0.7004900000000003, -0.6255299999999995, -0.5496499999999997, -0.4728999999999992, -0.3953799999999994, -0.31716999999999906, -0.2384000000000004, -0.1591799999999992,-0.11827, 0.0, 0.11827, 0.1591799999999992, 0.2384000000000004, 0.31716999999999906, 0.3953799999999994, 0.4728999999999992, 0.5496499999999997, 0.6255299999999995, 0.7004900000000003, 0.7744599999999995, 0.8474299999999992, 0.9193599999999993, 0.990219999999999, 1.060039999999999, 1.1288099999999996, 1.1965400000000006]
import matplotlib.pyplot as plt
import numpy as np

E = [1e3*i for i in E]
P = [i*0.248945227832799346 for i in P]

plt.figure()
plt.plot(P,E,'b.-')
plt.xlabel('Polarization ($C/m^2$)')
plt.ylabel('Energy ($meV$)')

f1 = np.polyfit(P, E, 6)
poly = ''
for i in range(len(f1)):
poly += '{}*x^{} + '.format(f1[i],len(f1)-i-1)
print(poly)

x = [ -0.35+0.7*i/99 for i in range(100)]
Eval = np.polyval(f1,x)
plt.figure()
plt.plot(x,Eval,'b')
plt.plot(P,E,'r*')
plt.xlabel('Polarization ($C/m^2$)')
plt.ylabel('Energy ($meV$)')
plt.show()

拉曼计算

VASP中计算拉曼谱主要使用vasp_raman.py这个脚本,见github:
https://github.com/raman-sc/VASP

使用raman-sc有几个注意的地方:https://zhuanlan.zhihu.com/p/354651066

作者这里提供了一个算例。供大家参考链接:https://pan.baidu.com/s/1dTCP_0bRcRo0uIy0rBGttQ
提取码:4vp8

  1. 这个代码已经年代久远,使用python2写的,对大部分新学vasp的同学来说可能装的anaconda3,既使用的是python3,与python2代码是不兼容的,作者这里对vasp_raman.py进行了改写,方便python3环境运行这段代码。下载见附件。
  2. raman.sub既提交代码的脚本,需要根据自己的软件环境进行更改,我这里给出我的提交脚本供大家参考:
    export VASP_RAMAN_RUN=’mpirun -np 64 /public/home/wangxiaohui/vasp/vasp.5.4.4/bin/vasp_std &> job.out ‘意思是运行vasp的命令,找一下vasp_std所在路径

python /public/home/wangxiaohui/WORKSPACE/Ti3C2Li2/raman/zhenraman/vasp_raman.py > vasp_raman.out这里是运行raman.py代码的意思,既raman.py代码所在路径

  1. 查看raman.py代码内容得知,POSCAR.phon和OUTCAR.phon是前一步计算的结果文件,不要被作者的例子给坑了。并且POSCAR.phon,OUTCAR.phon和INCAR等文件放在同一个文件夹下。

按照脚本VASP需要计算两步:

第一步,计算Gamma点声子。用有限位移方法或者DFPT方法。

第一步结束后,执行,

cp POSCAR POSCAR.phon

cp OUTCAR OUTCAR.phon

第二步,计算宏观介电常数的导数。用vasp_raman.py 脚本。

DFPT: LEPSILON=.TRUE.

频率依赖的介电矩阵: LOPTICS=.TRUE.

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#!/bin/bash
#PBS -l select=1:ncpus=32:mpiprocs=32
#PBS -l walltime=01:00:00
#PBS -q debug
#PBS -j oe
#PBS -N Example
#PBS -V

cd $PBS_O_WORKDIR

ulimit -s unlimited # remove limit on stack size

export VASP_RAMAN_RUN='mpirun -n 12 /u/afonari/vasp.5.3.2/vasp.5.3/vasp &> job.out' #a按照自己路径更改
export VASP_RAMAN_PARAMS='01_10_2_0.01' #b按照自己材料更改

python27 vasp_raman.py > vasp_raman.out

a. VASP_RAMAN_RUN设置调用VASP的命令

b. VASP_RAMAN_PARAMS=01_10_2_0.01

1
2
3
4
FIRST_MODE - 整数, 第一个计算的振动模
LAST-MODE - 整数, 最后一个计算的振动模
NDERIV - 整数, 有限位移的策略选择,目前只有2这一种选择
STEPSIZE - 浮点数, 有限位移步长,单位是Angstroms,一般选默认就可以

注意:以上这种设置格式允许我们并行计算大体系拉曼,我们可以在计算大体系的时候将Γ点的振动频率分为很多段分别计算,这样会大大提高我们计算大体系的效率。

在多段计算的时候。具体细节上要把每一段计算的OUTCAR.000n.-1.out和OUTCAR.000n.-1.out文件拷贝到同一个文件夹下,然后注释掉vasp_raman.py脚本中的运行VASP的命令。运行python2.7 vasp_raman.py > vasp_raman.out 即可得到整段的数据。

vasp_raman.py文件中执行VASP的行为331行:

1
2
3
else: # run VASP here
print "[__main__]: Running VASP..."
os.system(VASP_RAMAN_RUN) #注释掉这一行

加展宽

使用脚本broaden.py给出展宽谱,可以选择洛伦兹线型或者高斯线型。

$ python2 broaden.py result.txt
生成result.txt-broaden.dat文件,画图即可。

需要注意的是,broaden.py的第23和24行分别代表频率和强度所在列。其设置文件中对应的列。

cm1 = hw[:,1]
int1 = hw[:,2]

红外谱计算

  1. 红外谱主要是因为振动模在振动前后偶极矩发生变化,从而产生后外吸收。在振动过程中引起偶极矩改变的振动都具有红外活性。

  2. 根据上述原理,我们首先需要在计算过程中计算Gamma点的振动模,其次,我们需要计算伯恩有效电荷。

有限位移方法INCAR关键参数:

IBRION = 5
NWRITE = 3
POTIM = 0.015
LEPSILON = .TRUE
DFPT方法INCAR关键参数:

IBRION = 7
NWRITE = 3
NSW = 1
LEPSILON = .TRUE.
NELMIN=6

另外,建议进行k点测试直到伯恩有效常数收敛。

  1. 后处理

    计算得到OUTCAR,然后文件夹下运行脚本。运行结果如下:

sh IR.sh

上面的振动模和活性是归一化以后的结果,也就是运行结果intensities/results/results.txt的内容。文件intensities/results/exact.res.txt是归一化之前的内容。

  1. 展宽

使用脚本broaden.py给出展宽谱,可以选择洛伦兹线型或者高斯线型。

$ python2 broaden.py result.txt
生成result.txt-broaden.dat文件,画图即可。