10.1 铁电极化计算

发表于 更新于 分类于

介绍

铁电极化强度（ferroelectric polarization）是指在铁电材料中，由外加电场引起的电偶极矩的强度。铁电材料是一类具有自发极化特性的材料，即在没有外加电场时，这些材料内部的电偶极矩已经排列成某种方向。铁电极化强度描述了这种极化的强弱，它是材料的重要电学特性之一。

四种常见的铁电起源

  1. 离子偏移
    这种是最常见的铁电类型，比如BaTiO3中的Ti离子偏移中心位置；CuInP2S6中Cu离子上(下)偏移中心位置诱导的极化上(下)

  2. 极性分子团(polar molecular groups)比如有机-无机杂化材料的中一维（1D）和二维（2D）结构的杂化铁电材料，由极性分子诱导铁电。

  3. 电荷重布居(Charge redistribution)通过相邻层的占据态与邻近层的未占据态之间的杂化，可以实现层间电荷的重新分布，从而诱导出平面外的电偶极子，比如滑移铁电。

  4. 自旋比如正交晶系的TbMnO3，它具有非共线自旋结构的反铁磁性，通过逆Dzyaloshinskii Moriya相互作用产生净电极化。

下面介绍两种VASP计算铁电极化强度的方法。第一种，也是目前文献中最常见的方法，Berry phase方法。此外，还可以通过其他方法来计算。偶极矫正法是一种通过直接计算材料中偶极矩来估算铁电材料极化强度的简化方法。它通常用于快速评估铁电材料的极化强度，特别是在不需要高精度结果或没有使用现代极化理论（如 Berry 相位方法）的情况下。

例子1：离子型铁电

本文主要计算经典的铁电材料 BaTiO3 的极化强度。

原文链接:https://www.cnblogs.com/ghzhan/articles/16305679.html

BaTiO3 是钙钛矿结构，它的铁电相结构和中心对称结构如图所示，属于四方晶系。
该铁电材料是沿着 c 方向极化，主要是 O 离子和 Ti 离子的移动的贡献。 例如，上图中的左边铁电相极化方向朝下，右边铁电相极化方向朝上。极化强度详细计算步骤如下:

1 结构优化

先构造 BaTiO3 两种极化方向的晶格结构，并用 VASP 进行结构优化得到 CONTCAR；
铁电相POSCAR

1
2
3
4
5
6
7
8
9
10
11
12
13
Ti   Ba   O    
1.00000000000000       
 3.9944999999999999   0.0000000000000000   0.0000000000000000   
 0.0000000000000000   3.9944999999999999  -0.0000000000000000   
 0.0000000000000000   0.0000000000000000   4.0335000000000001
 Ti   Ba   O 
 1  1  3  
 Direct  
 0.5000000000000000  0.5000000000000000  0.5142700000000000    
 0.0000000000000000 -0.0000000000000000  0.0000000000000000    
 0.5000000000000000  0.5000000000000000  0.9744770000000000    
 0.5000000000000000  0.0000000000000000  0.4876180000000000    
 0.0000000000000000  0.5000000000000000  0.4876180000000000

1
2
3
4
5
6
7
8
9
10
11
12
13
Ti   Ba   O    
1.00000000000000        
3.9944999999999999   0.0000000000000000   0.0000000000000000   
0.0000000000000000   3.9944999999999999  -0.0000000000000000   
0.0000000000000000   0.0000000000000000   4.0335000000000001
Ti   Ba   O 
1  1  3  
Direct  
0.5000000000000000  0.5000000000000000  0.4857300000000000    
0.0000000000000000 -0.0000000000000000  0.0000000000000000    
0.5000000000000000  0.5000000000000000  0.0255230000000000    
0.5000000000000000  0.0000000000000000  0.5123820000000000    
0.0000000000000000  0.5000000000000000  0.5123820000000000

中心对称相POSCAR

1
2
3
4
5
6
7
8
9
10
11
12
13
Ti   Ba   O    
1.00000000000000        
3.9944999999999999   0.0000000000000000   0.0000000000000000   
0.0000000000000000   3.9944999999999999  -0.0000000000000000   
0.0000000000000000   0.0000000000000000   4.0335000000000001
Ti   Ba   O 
1  1  3  
Direct  
0.5000000000000000  0.5000000000000000  0.5000000000000000    
0.0000000000000000  0.0000000000000000  0.0000000000000000    
0.5000000000000000  0.5000000000000000  0.0000000000000000    
0.5000000000000000  0.0000000000000000  0.5000000000000000    
0.0000000000000000  0.5000000000000000  0.5000000000000000

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

INCAR要注意：
铁电极化计算建议设置高精度；DIPOL = 不要设置在原子及迁移路径上，设置在真空层的一边/质心。

批量化计算，计算程序根据相应计算环境调整。

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 方向极化的。对于具体的情况需要自行修改。离子与电子相加即可，然后用NP相-铁电相即可得到极化强度P。

注意

在铁电极化计算过程中经常会出现参考相为非半导体的情况，这种情况下可以：①镜像法：以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。

4 处理和分析数据。

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

1
2
1e =  1.602176634 * 10^-19 C  
1Å = 10^-10 m

三维体系的极化强度： 极化值除以体积。单位为 $e/Å^2$；

二维体系的极化强度：极化值除以面积。单位为 $e/Å$。

在这个例子中，BaTiO3 的体积为 64.3586 $Å^3$。

不同 image  的极化强度不是连续的，这是和选取的原胞有关，需要考虑极化量子的影响。BaTiO3 沿着 c 方向极化，所以需要对该极化值加减整数倍的 c 方向的晶格常数。

选取最靠近中间极化强度为 0 的那条曲线，即为 BaTiO3 的极化强度曲线：除以体积并进行简单的单位换算后为:另外，我们也可以绘制极化值 P 与能量 E 曲线

因此， BaTiO3 的极化强度大约为 0.248 $C/m^2$，与文献中的实验结果 0.26 $C/m^2$ 吻合。（Physical Review, 99(4), 1161–1165, 1955 ）

拟合 Landau-Ginzburg 公式
$$
E = \sum_{i} \left( \frac{A}{2} P_i^2 + \frac{B}{4} P_i^4 + \frac{C}{6} P_i^6 \right) + \frac{D}{2} \sum_{i,j} (P_i - P_j)^2
$$

分别计算各个原子位移（原子移动方式可以类似于线性插值法）后结构的极化值，通过L-G公式拟合极化与能量得到前三项的系数A,B和C。由具有铁电性的原胞（FE结构）建立441超胞，计算超胞能量E1。然后将441超胞中其中一个原胞的FE结构换成-FE结构（-FE结构），代表一个dipole发生了翻转，计算出超胞能量E2，根据朗道有效哈密顿量在这两种情况下的表达式，并比较能量E1和E2，就可以得到D的数值。

拟合该多项式曲线，得到$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
$$

 绘制能量曲线的脚本

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曲线，以及考虑极化量子的脚本:

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 曲线的脚本:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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()

参考:

https://chengcheng-xiao.github.io/post/2019/12/25/Wannier_center_polarization.html

https://www.bilibili.com/read/cv14756243/

https://zhuanlan.zhihu.com/p/358517335

https://zhuanlan.zhihu.com/p/537595648

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

https://zhuanlan.zhihu.com/p/445884561

5). 居里温度计算的两种方法

  1. AIMD方法计算居里温度Tc（计算量巨大，准确）使用AIMD分别计算不同温度下的超胞，通过统计学得到平均极化距离d，通过d找到对应温度下的极化值，然后通过slogical1函数拟合得到$T_c$。
  2. MC方法计算居里温度$T_c$（计算速度极快，建议）通过L-G公式进行拟合即可，代码及教程详细可看https://github.com/Chengcheng-Xiao/mpiPyMC

滑移铁电

华中科技大学的吴梦昊教授，在2017年首次提出层间滑移铁电，通过两层之间的滑移，实现电极化翻转。

下面以双层GdI2为例，计算这个材料的铁电极化强度。

Xun W, Wu C, Sun H, et al. Coexisting magnetism, ferroelectric, and ferrovalley multiferroic in stacking-dependent two-dimensional materials[J]. Nano letters, 2024, 24(11): 3541-3547.

1 首先准备优化好的铁电相，即AB堆叠的POSCAR：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
CONTCAR                                    
1.00000000000000          
4.1426458804503365    0.0000000000000000    0.0000000508502232    
-2.0713229953957617    3.5876365788074991   -0.0000000399383119     
0.0000004448876875    0.0000000102694673   30.0000000000000000   
I    Gd     
4     2
Direct  
0.3333333244151952  0.6666666877589731  0.2311447203064440  
0.3333333924379757  0.6666666704222977  0.0889274910801614  
0.6666666411760427  0.3333333528368931  0.4864429969506502  
0.6666666937613037  0.3333333375676722  0.3442418595816502  
0.6666666947651289  0.3333333427541252  0.1601182205842356  
0.0000000004443691  0.0000000116600432  0.4152546984968589

在INCAR中加入参数：

1
2
DIPOL = 0.5 0.5 0.3
LCALCPOL = .TRUE.   #计算铁电极化开关

2 自洽计算

计算完成后，在OUTCAR中查找关键字dipole：由于我们关注的是c方向，因此只要看最后一列，用-397.00967 + 7.0125 = -389.99717，即为铁电相的极化值。

3 接下来计算顺电相，即AA堆叠，

优化好的POSCAR如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
CONTCAR                                    
1.00000000000000          
4.1426458804503365    0.0000000000000000    0.0000000508502232    
-2.0713229953957617    3.5876365788074991   -0.0000000399383119     
0.0000004448876875    0.0000000102694673   30.0000000000000000   
I    Gd     
4     2
Direct  
0.6756282937225784  0.3441349445403583  0.2202487956955387  
0.6756257399639071  0.3441404853191938  0.0779182529604372  
0.6708730704449056  0.3402390499976288  0.4974588775154418  
0.6708673133601825  0.3402448262892301  0.3551273278370464  
0.0089448210006132  0.0107926773902162  0.1491235435481827  
0.0042207755078186  0.0069279784633755  0.4262532274433494

同样在自洽的INCAR中加入参数：

1
2
DIPOL = 0.5 0.5 0.3
LCALCPOL = .TRUE.  #计算铁电极化开关

计算完成后查看顺电相的极化：

极化值为：-397.00991 + 7.00991 = -390

4 数据处理

铁电相减去顺电相，得到的值为：-389.99717-(-390)=0.00283，单位是电子/埃。

极化强度是单位面积上的极化值，进行单位转换，最终结果3.0467x10-13 C/m。

对比文献中的结果：数值还是非常接近的，差距可能在结构优化，或者参数设置的一些细节方面。至于为什么文献中是-12次方，我的是-13次方，我特地联系了通讯作者询问了一下，是他们笔误搞错了，正确结果就是-13次方。最后可能有人发现了，似乎没用到极化量子，由于计算的结果比较凑巧，刚好在同一个极化量子的范围内。一般情况下，可能算出来的极化值会相差很多个极化量子。

原链接：https://mp.weixin.qq.com/s/1j97g5A1Ny95EOn3cqKh-A

功函数和静电势计算

VASP 中可使用 LVHAR 参数进行功函数计算控制，输出文件中 LOCPOT 就是我们想要的静电势贡献，格式和 CHGCAR 一样，需计算 Z 方向的平均贡献,利用vaspkit -task 422/426.

VASP 中直接使用 LDIPOL 和 IDPOL 即可开启它的偶极校正功能

1
2
3
4
LDIPOL = .TRUE. 表示打开偶极校正
IDIPOL = 3 表示偶极校正方向为第三晶格矢
DIPOL = 0.5 0.5 0.5 表示体系的中心，以分数坐标表示
LVHAR = .TRUE.

真空能级在打开 LVHAR 时可能去读取 OUTCAR 直接得到

grep vacuum OUTCAR

1
# vacuum level on the upper side and lower side of the slab         2.807         3.188

注意 OUTCAR 中的真空能级需进行费米能级修正，即减去 OUTCAR 中的 E-fremi

注意

功函数：将一个固体内部的电子移动到真空所需的能量。（类似于逸出功）。

真空能级：固体表面外真空中自由电子所具有的能量。换句话说，电子跑出固体表面并达到这个能级后即可认为它自由。

  1. 因 VASP 所适用的体系是周期性体系，使用它来模拟实验中的 Slab 模型时会取一个相当大的真空层来隔绝相邻两个周期中 Slab 的相互作用。理想情况下，真空层中的功函数应当是一条水平的直线（函数值为定值）。但如果表面的两侧并非对称，即其中一侧吸附了分子时， 这两侧的功函数存在差异，此时如果不进行偶极校正，真空中的功函数会是一条斜线；而经过偶极校正后，功函数会出现一个阶梯，阶梯两侧附近的曲接近水平。

  2. 在 VASP 中 LVHAR 参数可以使 VASP 输出体系的功函数文件 LOCPOT 。LOCPOT 文件本身是 Volumetric data ，它的格式与 CHGCAR 一样。一般而言，用户关心的功函数是垂直于表面方向上的数据，因此在得到 LOCPOT 后需要对 xy平面内的数据做平均，然后乘以晶胞的体积，就得到我们需要的功函数信息。

  3. 若体系难以收敛可手动计算原子坐标平均值并设置 DIPOL = 0.5 0.5 ，若修正效果不到预期可以打开偶极校正并设置上述 DIPOL = 0.5 0.5 进行弛豫计算；
    关闭 DIPOL，再次进行弛豫计算；
    关闭 DIPOL，打开 LVHAR = .TRUE. 进去静态计算得到功函数；
    若出现功函数台阶即为计算成功。

  4. 如果计算失败请检查 POSCAR 中晶格基矢对应分量是否为 0。如计算第三个方向偶极修正请保证前两个基矢的 z 分量为零且第三个基矢只有 z 分量不为零。另外要注销并行参数。