4.2 光学性质
介电常数
首先进行结构优化,参考官网: https://www.vasp.at/wiki/index.php/Dielectric_properties_of_SiC
静态介电常数:
方法1:DFPT方法;
1
2
3
4
5IBRION = 8
NSW=1
LEPSILON = .TRUE.
LRPA=.TRUE.#RPA考虑局域场效应,默认关闭,即IP近似
LPEAD=.TRUE.#有限差分法
方法2:Response to finite electric fields
1 | LCALCEPS = .TRUE. |
在OUTCAR中抓取:
1 | ```MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))``` |
frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects)
andfrequency dependent REAL DIELECTRIC FUNCTION (independent particle, no local field effects)
Note:对于光学性质的计算,也就是计算材料的介电函数,需要足够多的空带和致密的K网格点,使其达到非常好的收敛状态,我们才可以得到合理的光学性质;因此通常计算中,一般设置NBANDS为默认值(可在自洽的OUTCAR中以关键字NBANDS查找到)的2-3倍,K网格点为自洽值或适当增加。
使用optics.sh脚本得到介电函数的实部和虚部
1 | # extract image and real parts of dielectric function from vasprun.xml |
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
极化强度
- 极化强度计算
首先选取极化相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
20SYSTEM = 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. - 收集结果: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 是钙钛矿结构,它的铁电相结构和中心对称结构如图所示,属于四方晶系。
先构造 BaTiO3 两种极化方向的晶格结构,并用 VASP 进行结构优化得到 CONTCAR;
将上一步优化后的两个结构分别放入创建好的 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 文件。
对每个文件夹的结构进行一次 VASP 自洽运算,INCAR 文件里面需要额外设置 DIPOL 和 LCALPOL 参数。
DIPOL 参数可以选取任一坐标,但是保证同一体系采用相同值。
LCALPOL 参数是打开 极化运算。
INCAR 文件:
1
2
3
4
5
6
7
8
9
10
11
12SYSTEM=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计算完成后运行命令:
1
2
3grep '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 方向极化的。对于具体的情况需要自行修改。
处理和分析数据。
首先要理清数据单位。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) 的总能:
用dipole_c.dat 文件 (单位: e*Å):不同 image 的极化值如图所示
可以发现,该极化强度不是连续的,这是和选取的原胞有关,需要考虑极化量子的影响。BaTiO3 沿着 c 方向极化,所以需要对该极化值加减整数倍的 c 方向的晶格常数。这样我们得到下图:选取最靠近中间极化强度为 0 的那条曲线,即为 BaTiO3 的极化强度曲线:
除以体积并进行简单的单位换算后为:
另外,我们也可以绘制极化值 P 与能量 E 曲线
因此, 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$
- 绘制能量曲线的脚本
{.line-numbers} 1
2
3
4
5
6
7
8
9
10
11with 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曲线,以及考虑极化量子的脚本:
1 | with open('dipole_c.dat','r') as f: |
绘制极化值 P 与能量 E 曲线的脚本:
1 | 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] |
拉曼计算
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
- 这个代码已经年代久远,使用python2写的,对大部分新学vasp的同学来说可能装的anaconda3,既使用的是python3,与python2代码是不兼容的,作者这里对vasp_raman.py进行了改写,方便python3环境运行这段代码。下载见附件。
- 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代码所在路径
- 查看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.
1 | #!/bin/bash |
a. VASP_RAMAN_RUN设置调用VASP的命令
b. VASP_RAMAN_PARAMS=01_10_2_0.01
1 | FIRST_MODE - 整数, 第一个计算的振动模 |
注意:以上这种设置格式允许我们并行计算大体系拉曼,我们可以在计算大体系的时候将Γ点的振动频率分为很多段分别计算,这样会大大提高我们计算大体系的效率。
在多段计算的时候。具体细节上要把每一段计算的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 | else: # run VASP here |
加展宽
使用脚本broaden.py给出展宽谱,可以选择洛伦兹线型或者高斯线型。
$ python2 broaden.py result.txt
生成result.txt-broaden.dat文件,画图即可。
需要注意的是,broaden.py的第23和24行分别代表频率和强度所在列。其设置文件中对应的列。
cm1 = hw[:,1]
int1 = hw[:,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点测试直到伯恩有效常数收敛。
后处理
计算得到OUTCAR,然后文件夹下运行脚本。运行结果如下:
sh IR.sh
上面的振动模和活性是归一化以后的结果,也就是运行结果intensities/results/results.txt的内容。文件intensities/results/exact.res.txt是归一化之前的内容。
- 展宽
使用脚本broaden.py给出展宽谱,可以选择洛伦兹线型或者高斯线型。
$ python2 broaden.py result.txt
生成result.txt-broaden.dat文件,画图即可。