4.2 光学性质

发表于 2023-05-11 更新于 2023-12-26 分类于 DFT学习

介电常数

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

静态介电常数：

方法1：DFPT方法;

{.line-numbers}

IBRION = 8
NSW=1
LEPSILON = .TRUE.
LRPA=.TRUE.#RPA考虑局域场效应，默认关闭，即IP近似
LPEAD=.TRUE.#有限差分法

方法2：Response to finite electric fields

{.line-numbers}

LCALCEPS = .TRUE.
IBRION=6
NFREE=2
ISMEAR =  0
SIGMA  =  0.01

在OUTCAR中抓取：

STATIC DIELECTRIC TENSOR (including local field effects in DFT)``` 电子贡献或者

```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}

# 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

极化强度

极化强度计算
首先选取极化相FE（PhaseB, P, 铁电相）和参考相NP（PhaseA, P=0, 中心对称相），分别以下列INCAR计算，

{.line-numbers}

 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.

收集结果：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 ，属于四方晶系。

先构造 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
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
计算完成后运行命令：
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 方向极化的。对于具体的情况需要自行修改。
处理和分析数据。

首先要理清数据单位。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$
0408

绘制能量曲线的脚本

{.line-numbers}

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}

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}

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

这个代码已经年代久远，使用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.

{.line-numbers}

#!/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

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]

红外谱计算

红外谱主要是因为振动模在振动前后偶极矩发生变化，从而产生后外吸收。在振动过程中引起偶极矩改变的振动都具有红外活性。
根据上述原理，我们首先需要在计算过程中计算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文件，画图即可。