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