3.1 热输运1 声子谱

发表于 更新于 分类于

声子谱

声子谱的计算主要有两种方法，一种是直接法，另一种是微扰密度泛函方法（DFPT）。

  • 直接法，或称frozen-phonon方法，是通过在优化后的平衡结构中引入原子位移，计算作用在原子上的Hellmann-Feynman力，进而由动力学矩阵算出声子色散曲线。直接法的缺陷在于它要求声子波矢与原胞边界supersize
    正交，或者原胞足够大使得Hellmann-Feynman力在原胞外可以忽略不计。这使得对于复杂系统，如对称性高的晶体、合金、超晶格等材料需要采用超原胞。超原胞的采用使计算量急剧增加，极大的限制了该方法的使用。
  • 1987年，Baroni、Giannozzi和Testa提出了DFPT方法。DFPT通过计算系统能量对外场微扰的响应来求出晶格动力学性质。该方法最大的优势在于它不限定微扰的波矢与原胞边界正交，不需要超原胞也可以对任意波矢求解。Castep、Vasp等采用的是一种linear response theory 的方法（或者称为 density perturbation functional theory，DFPT），直接计算出原子的移动而导致的势场变化，再进一步构造出动力学矩阵。进而计算出声子谱。

DFPT计算步骤

  • 结构优化
    声子谱的计算需要对原胞结构做高精度的充分优化，否则很容易出现虚频；
    {.line-numbers}
    1
    2
    EDIFF  = 1E-08
    EDIFFG = -1E-06
    可以采用分步优化，直到达到精度。然后优化后的结果，cp CONTCAR POSCAR：

$ phonopy -d --dim="4 4 1"

之后会生成一堆文件，其中有 SPOSCAR 和POSCAR-00* ，分别对应不同的方法密度泛函微扰法 (DFPT) 以及有限位移法。

{.line-numbers}
1
2
3
4
5
6
$ ls
phonopy_disp.yaml  POSCAR-002  POSCAR-005  POSCAR-008  POSCAR-011  POSCAR-014  SPOSCAR
POSCAR             POSCAR-003  POSCAR-006  POSCAR-009  POSCAR-012  POSCAR-015
POSCAR-001         POSCAR-004  POSCAR-007  POSCAR-010  POSCAR-013  POSCAR-016
$ cp POSCAR POSCAR-unitcell
$ cp SPOSCAR POSCAR

KPOINTS需适当减小，可以的话最好再进行一次收敛测试
INCAR
{.line-numbers}
1
2
3
4
5
6
7
8
9
10
IBRION =  8            (determines the Hessian matrix using DFPT)
EDIFF  =  1E-08        (SCF energy convergence; in eV) 
PREC   =  Accurate     (Precision level)  
ENCUT  =  500          (Cut-off energy for plane wave basis set, in eV) 
IALGO  =  38           (Davidson block iteration scheme)
LREAL  = .FALSE.       (Projection operators: false)
LWAVE  = .FLASE.       (Write WAVECAR or not)
LCHARG = .FLASE.       (Write CHGCAR or not) 
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence) 
NSW    = 1

上面参数 IBRION = 8 NSW = 1是注意修改的。
提交VASP计算

后处理
phonopy --fc vasprun.xml

读取vasprun.xml，生成FORCE_CONSTANTS文件。

编辑band.conf文件，该文件给出了高对称点路径的信息
Note:

{.line-numbers}
1
2
3
4
5
ATOM_NAME = * *  #元素名称，顺序与POSCAR保持一致
DIM = 4 4 1      #建立超胞的大小，--dim后面的
BAND = 0.5 0.0 0.0  0.0 0.0 0.0  0.333333 0.333333 0.0  0.5 0.0 0.0 #高对称点路径，每组高对称点之间用两个空格隔开
FORCE_CONSTANTS = READ #表示读取力常数文件FORCE_CONSTANTS
FC_SYMMETRY = .TRUE.

准备mesh.conf文件，如下所示：
{.line-numbers}
1
2
3
ATOM_NAME = Si
DIM = 2 2 2
MP = 24 24 24

{.line-numbers}
1
2
3
4
phonopy --dim="4 4 1" -c POSCAR-unitcell -p -s band.conf
(单位是THz)或者
phonopy --dim="4 4 1" --factor=524.147 -c POSCAR-unitcell -p -s band.conf
(单位是cm-1)

（读取band.conf、FORCE_CONSTANTS、POSCAR-unitcell，输出band.yaml、phonopy.yaml以及pdf图）

1
2
3
phonopy-bandplot  --gnuplot> phono.dat
如果上述命令不行，可以用旧版本phonopy 
bandplot --gnuplot> phonon.dat

读取band.yaml，就可以获得二阶力常数矩阵声子谱的数据，其中第一行为高对称点的坐标。

准备mesh.conf以计算态密度和PDOS：

{.line-numbers}
1
2
3
4
5
6
ATOM_NAME = 
DIM = 4 4 1
MP=12 12 1
PDOS=1 2, 3 4, 5 6
FORCE_CONSTANTS = READ
BAND_POINTS = 71

运行命令：phonopy –dim='4 4 1' -c POSCAR-unitcell mesh.conf或者乘上因子修改单位。phonopy –dim='4 4 1' --factor=524.147 -c POSCAR-unitcell mesh.conf
得到projected_dos.dat文件，拖到Origin中画图，得到声子态密度图像（从第二列数据开始，每一列对应每个原子的态密度），对比着之前计算得到的声子谱可以分析声学支与光学支之间的相互作用关系

有限位移方法

这里需要准备和 POSCAR-00* 数目一样的文件夹，并把 POSCAR-00* 复制到每个文件夹中并命名为 POSCAR 。每个文件夹内都需要准备 POTCAR 、 POSCAR 、 INCAR 、 KPOINTS 。并且分别计算。
可用以下sh脚本快速进行：手动修改第一行文件夹数目和vasp运行文件名称。

{.line-numbers}
1
2
3
4
5
6
7
8
9
for i in {001..006}
do 
mkdir $i
cd $i
cp ../POSCAR-$i POSCAR
cp ../{INCAR,KPOINTS,POTCAR,sub.sh} .
运行提交任务命令
cd ../
done

INCAR 可以设置为如下内容：
INCAR设置如下(静态计算)：

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
PREC = Accurate
IBRION = -1
NSW = 1
ENCUT = 500
EDIFF = 1.0e-08
EDIFFG = -0.001
ISMEAR = 0
SIGMA = 0.05
IALGO = 38
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.

里面IBRION = -1 NSW = 1注意修改

KPOINTS需适当减小，可以的话最好再进行一次收敛测试
提交VASP计算

后处理

{.line-numbers}
1
phonopy -f ./disp-*/vasprun.xml

准备band.conf

{.line-numbers}
1
2
3
4
5
ATOM_NAME = In Se
DIM = 4 4 1
BAND = 0.5 0.0 0.0  0.0 0.0 0.0  0.333333 0.333333 0.0  0.5 0.0 0.0
FORCE_CONSTANTS = WRITE  #注意这里是 WRITE
FC_SYMMETRY = .TRUE.
{.line-numbers}
1
2
3
4
5
6
7
8
phonopy -c POSCAR-unitcell band.conf -p -s
phonopy --factor=524.147 -c POSCAR-unitcell -p -s band.conf
#旧版本phonopy
bandplot  --gnuplot> phonon.dat
就可以获得二阶力常数矩阵 FORCE_CONSTANTS 和声子谱的数据
#新版本phonopy
phonopy-bandplot --gnuplot > phonon.dat

Mind

  • 如果计算完的声子谱其他地方都好，就是在Γ点有一点点虚频，那么这个材料很可能是稳定的，只是你优化做得不够好，进一步提高优化的精度可消除这一点点虚频。
  • 对于二维材料，如果在Γ点出现很小的虚频，基本可以认为这个材料是稳定的，大部分二维材料都会有此现象；尤其是VASP结合phonopy code计算二维材料的声子谱在Γ点更是容易出现虚频；使用Quantum-ESPRESSO的PWSCF和PH模块计算声子谱对于内存的需求较小一些，且对于二维材料的声子谱计算更友好一些，后期会详细介绍Quantum-ESPRESSO计算声子谱的具体步骤。
  • 声子谱是表示组成材料原子的集体振动模式。如果材料的原胞包含n个原子，那么声子谱总共有3n支，其中有3条声学支，3n-3条光学支。声学支表示原胞的整体振动，光学支表示原胞内原子间的相对振动。

计算出的声子谱有虚频，往往表示该材料不稳定。 有些情况下，我们可以利用虚频信息使不稳定的材料变得稳定。如果声子谱的一条声学支存在虚频，主要位于Γ点和M点1/2处（对应倒格矢的1/4位置）。倒格矢的1/4，对应晶格长度的4倍。我们可能需要将原胞沿上述倒格矢方向扩大四倍，进一步优化原子位置，才可能得到比较稳定的晶胞。

消虚频

  1. DFPT方法与有限位移法互换尝试。或者更换标准胞扩胞计算。
  2. 扩胞，在虚频很小时可以尝试扩大扩胞倍数。
  3. 提高精度优化后计算声子谱。声子谱对力的精度要求特别高，所以提高精度多次优化后再计算声子谱也是可行的。
  4. 修改POTCAR，尝试不同的赝势方法。（不太建议，因为需改之后所有计算都需要统一成稳定构型的赝势。不过有效）
  5. 修改POSCAR原子坐标，再优化，再计算声子谱。
    修改方法：对于DFPT计算的声子谱，使用grep cm-1 OUTCAR查看频率，寻找最大的虚频频率。然后在OUTCAR中从最后往前找到最大的虚频频率对应的坐标。
    dx dy dz对应虚频的位移，把这些坐标与POSCAR中的坐标相加或者0.1-0.9之间的倍数，形成新的POSCAR。进行优化，再计算声子谱。
  6. 不排除真的算不出来没虚频的，据我所知，SnTe就没见到谁大大方方发文章谈声子谱的（算出来虚频消除不了不好意思讲）。
  7. INCAR参数对虚频的影响：ADDGRID影响很小，LREAL有影响。截断能和K点密度在收敛情况下影响小。

结构优化消虚频

X, Y, Z为体系中原子的坐标，也就是POSCAR中的内容。 dx, dy,dz为虚频对应的原子振动方向上的位移。我们将（dx, dy,dz）这个位移乘以一个介于0-1之间的校正因子，然后跟POSCAR中的坐标加起来即可。以z方向为例，校正因子设为0.1， 微调后的z方向坐标为： 10.709980 + （-0.676634）0.1。
脚本vib_correct.py如下：https://www.bigbrosci.com/2022/04/29/A29/

{.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
40
41
42
43
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py 
This is a temporary script file.
"""
import numpy as np
from ase.io import read, write
import os

#获取虚频开始行
l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
    lines = f_in.readlines()
    wave_num = 0.0
    for num, line in enumerate(lines):
        if 'f/i' in line:
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num: #获取最大的虚频
                wave_num = wave_tem
                l_position = num+2
# ASE读POSCAR
model = read('POSCAR')
model_positions = model.get_positions()
num_atoms = len(model)
#print(model_positions)

# 获取虚频对应的OUTCAR部分 
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行
#print(vib_lines)
vib_dis = []
for line in vib_lines:
    #model_positions = [float(i) for i in line.rstrip().split()[:3]]
    vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列
    vib_dis.append(vib_infor)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。
model.positions = new_positions

###保存结构
write('POSCAR_new', model, vasp5=True) # POSCAR_new是微调后的结构，用于下一步的计算（别忘了把POSCAR_new改成POSCAR）。

过渡态虚频

如果过渡态计算的时候出现好几个虚频，自己直接写个python脚本也能实现。不过使用ASE可以很方便地获取POSCAR中原子数目，保存原子固定的信息。参考如下。

{.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
40
41
42
43
44
45
46
47
48
49
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
import os

#获取虚频开始行
l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
    lines = f_in.readlines()
    wave_num = 0.0
    for num, line in enumerate(lines):
        if 'f/i' in line:
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num: #获取最大的虚频
                wave_num = wave_tem
                l_position = num+2
# 读POSCAR
with open('POSCAR') as f_pos:
    lines_pos = f_pos.readlines()

# 获取虚频对应的OUTCAR部分    
num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()])
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行

model_positions = []
vib_dis = []
for line in vib_lines:
    position = [float(i) for i in line.rstrip().split()[:3]]
    vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列
    model_positions.append(position)
    vib_dis.append(vib_infor)

# 微调结构
model_positions = np.array(model_positions)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。

###保存结构
f_out = open('POSCAR_new','w')
f_out.writelines(lines_pos[:8])
f_out.write('Cartesian\n')
for i in new_positions:
    f_out.write(' '.join([str(coord) for coord in i]) + '  F  F  F\n')
f_out.close()

phonopy声子群速度

操作步骤

  1. 首先对结构高精度优化后扩胞，计算得到没有虚频的声子谱；
  2. 对band.conf文件进行修改，并添加群速度计算相关参数GROUP_VELOCITY = .TRUE.和GV_DELTA_Q = 0.01，程序需要读取本征值和本征矢，所以必须设置 EIGENVECTORS =.TRUE.，此外为了防止vaspkit在导出数据时报错，需要删除BAND_LABELS高对称路径名称，修改后band.conf文件内容如下；

而且具体使用过程中，关于执行命令还有执行文件的tag内容还有一些可替换的关系,–gv, –group_velocity 可在执行计算命令时替代执行文件中的GROUP_VELOCITY = .TRUE.内容， –gv_delta_q 可在执行计算命令时替代执行文件中的(GV_DELTA_Q)内容

  1. 对声子谱重新进行计算，打开phonopy.yaml文件可以看到，phonopy在计算声子谱的过程中考虑到了群速度；打开band.yaml文件中也保存了有关群速度的相关信息；同时在计算声子谱中phonopy所储存的数据文件mesh.yaml和band.yaml文件中也保存了有关群速度的相关信息。
  2. 随后进行数据提取，以vaspkit1.12版本为例，打开vaspkit，输入73，打开VASP2other Interface 子功能，选择739 Sort Phonon Band Structure for Phononpy，运行完成即可获得group velocity相关数据以及各原子的声子谱数据；在运行过程中我们可以看到软件说明find group velocity并写入输出文件中。
    1.5版本为789功能。

格林艾森常数（Grüneisen parameters）

为表示声子非简谐效应的参数。phonopy计算模式格林艾森参数方式见官网。
以NaCl为例，首先需要进行三个声子谱计算，分别是平衡体积（orig），稍小体积（minus）和稍大体积（plus）的结构，计算时改变晶格参数使得体积获得改变。

计算完成后获得不同体积下的原子间力常数文件FORCE_SETS 或 FORCE_CONSTANTS （后续命令中需要添加–readfc参数）。

随后通过phonopy-gruneisen（旧版gruneisen）命令获得对应于声子带结构中倒易空间中的路径的模式格林艾森常数（具体设定可根据自身计算情况更改）
phonopy-gruneisen orig plus minus --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --band="1/2 1/4 3/4 0 0 0 1/2 1/2 1/2 1/2 0.0 1/2" -p -c POSCAR-unitcell
在倒易网格中的模式格林艾森常数（具体设定可根据自身计算情况更改）：
phonopy-gruneisen orig plus minus --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --mesh="20 20 20" -p -c POSCAR-unitcell --color="RB"
phonopy-gruneisen 命令选项包括–dim，–mp, –mesh，–band，–pa, –primitive_axis，–readfc，–band_points，–nac，–factor，–nomeshsym，-p，-c，-s, –save，-o等，功能上均与phonopy命令相同，但除上述列出的选项外，可以在phonopy命令上使用的其它选项是不能在phonopy-gruneisen 上使用的。如有必要，可以自行修改phonopy-gruneisen代码。

如查看帮助，可输入
phonopy-gruneisen -h
通过 phonopy-gruneisenplot 命令导出计算得到的数据，使用方式类似于phonopy-bandplot。

原文链接：https://blog.csdn.net/icehoqion/article/details/130550807

热学性质

参考：
http://supersunsir.com/2020/08/01/vasp/vasp_phonopy/

热膨胀系数

我们通过准简谐近似（QHA）的方法计算热膨胀系数和等压热容等性质，

  • 首先，我们要把优化好的POSCAR加应变（比如0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04的缩放系数，这样就有9个结构文件），
  • 然后 把每一个POSCAR都去重复步骤（2）也就是去计算每一个缩放系数声子谱。大概就是这个效果
  • 最后我们通过phonopy里面的QHA模块 去得到此材料的非偕性质phonopy-qha -p -s e-v.dat *00/thermal_properties.yaml就可以得到很多热力学性质比如热膨胀系数，格林爱森参数，自由能，等压热容等等。
  1. 先进行高精度的晶格优化。
  2. 修改缩放系数为0.95-1.05新建10个POSCAR.
  3. 对11个POSCAR晶体进行声子谱计算
  4. thermo.conf文件后处理生成11个对应的thermal_properties.yaml
  5. QHA计算热膨胀系数

另一个写的步骤：

  1. 高精度的原胞优化
  2. 在此基础上施加形变（梯度可以设置为千分之一）
  3. 优化多个结构
  4. 使用DFPT方法计算每个结构的声子谱（如果结构有问题，使用phonopy –symmetry）
  5. 在每个文件夹中得到thermal_properties.yaml
  6. 运行脚本得到e-v.dat
  7. 运行phonopy-qha的相关命令得到各种dat文件
  8. Over

可以参考作者的thermo.conf文件：

{.line-numbers}
1
2
3
4
5
6
7
DIM = 2 2 2
ATOM_NAME = sio
MP = 31 31 31
TPROP = .TRUE.
TMAX = 2000
TSTEP = 2
FORCE_CONSTANTS = READ

然后是如何写e-v.dat文件

  • 大概意思就是找到之前计算的那11个POSCAR的vasprun.xml结果文件，然后记事本打开搜索最后一个e_wo_entrp。
  • 下载phonopy官方文件中Si-QHA的example，
  • 解压vasprun_xmls.tar.lzma我们可以得到算例中用到的vasprun.xml文件。

不同温度下热力学性质计算（自由能、热容、振动熵）

  • 在DFPT或者Finite displacement方法计算声子谱的基础上，得到mesh.yaml文件

  • 新建一个文件夹，将FORCE_CONSTANTS、POSCAR、mesh.conf放入该文件夹中（POSCAR为原胞）

  • 得到mesh.yaml文件

    {.line-numbers}
    1
    phonopy mesh.conf -c POSCAR

    mesh.conf的内容

    {.line-numbers}
    1
    2
    3
    4
    ATOM_NAME = Hf Si O
    DIM = 2 2 2
    MP = 13 13 13
    FORCE_CONSTANTS = READ

    得到thermal_properties.yaml文件

phonopy -c POSCAR mesh.conf -t
得到thermal_properties.dat文件
phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat
将dat文件导入到Origin中画图

dat文件的5列数据分别为：

{.line-numbers}
1
2
3
4
5
temperature[K]
Free energy [kJ/mol]
Entropy [J/K/mol]
Heat capacity [J/K/mol]
Energy [kJ/mol]

下面是能用到的各种命令

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
扩胞用
phonopy -d --dim="2 2 2" -c POSCAR_UNIT
查找晶格对称性
phonopy --symmetry
收集力常数
phonopy --fc vasprun.xml
得到band.yalm文件（好像不需要这一步）
phonopy --dim="2 2 2" -c POSCAR-UNIT -p band.conf
得到thermal_properties.yaml文件
phonopy -c POSCAR mesh.conf  -t

下面是产生e-v.dat文件的脚本，需要根据自己的目录名进行修改：

{.line-numbers}
1
2
3
4
5
6
7
8
9
for i in -0.005 -0.004 -0.003 -0.002 -0.001 0.001 0.002 0.003 0.004 0.005
do
 cd AC_Uni$i;
 pwd;
 volume=`grep volume vasprun.xml | tail -1| awk '{printf("%20.13f", $3)}'`;
 energy=`grep 'E0' OSZICAR | tail -1 | awk '{printf "%12.4f \t", $5}'`;
 echo ${volume} ${energy} >> ../e-v.dat;
 cd ../; #Exit the Strain folder
done

声子谱振动可视化

  1. 上面生成band.yaml文件，利用JMOL打开。https://jmol.sourceforge.net/
    或者网站：https://henriquemiranda.github.io/phononwebsite/phonon.html

  2. band.conf添加标签ANIME=0 0 0将000这个Q点的振动模式输出到anime.ascii,用v_sim软件打开即可。https://l_sim.gitlab.io/v_sim/index.en.html

  3. 运行命令 phonopy -f vasprun.xml-{001,002}得到FORSETS
    运行phonopy –dim=”3 3 3” –ct=”0 0 0”得到gamma点的振动模式密度

  4. DFPT方法得到的OUTCAR可以利用PyVibMS、JMOL

  5. 从band.yaml生成vesta读取：https://github.com/AdityaRoy-1996/Phonopy_VESTA

注意不要输出声子群速度，写出声子本征矢量。需要POSCAR.vesta格式。

从OUTCAR生成vesta读取：https://github.com/Stanford-MCTG/VASP-plot-modes

  1. 其他参考：https://www.bilibili.com/read/cv11848434/

声子谱反折叠

https://github.com/yuzie007/upho

热位移

热学性质是根据倒易空间采样网格上的声子频率计算得出的，必须设置Mesh tag进行使用，并且必须检查它们相对于取样网格的收敛性。 通常这种计算对计算量要求不高，因此随着采样网格密度的增加，收敛性很容易达到。此外，计算时Phonopy对频率的默认单位为Thz。
在DFPT和冷冻声子的计算基础上，在thermal_properties.conf 中设置：

1
2
3
4
5
6
TPROP = .TRUE. #进行热学性质计算
DIM = 1 1 1 #扩胞情况
MP = 16 16 8 #Mesh网格密度
TMAX = 2000
TMIN = 0
TSTEP = 10 #温度上下限默认为100,0，step默认值为10

计算得到thermal_properties.yaml文件：

phonopy -c POSCAR thermal_properties.conf -t

得到thermal_properties.dat文件：

phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat

五列数据分别为：temperature[K]、Free energy [kJ/mol]、Entropy [J/K/mol]、Heat capacity [J/K/mol]、Energy [kJ/mol]。

Phonopy官网还给出了两个参数：
PRETEND_REAL = .TRUE.
这里将虚频声子当成实声子进入计算，仅在有小虚频实作测试用；
CUTOFF_FREQUENCY = 0.1
默认值为0，将虚频声子频率视为负值，计算高于截止频率的声子。

热位移计算

计算投影到笛卡尔轴上的均方位移与温度的函数关系。声子频率以 THz 为单位（phononopy 的默认设置），用于获得均方位移。同样需要与Mesh tag共同使用，温度设置与前述相同。新建thermal_displacements.conf用于输入tag。

1
2
3
4
5
6
TDISP = .TRUE. #进行均方位移计算
DIM = 1 1 1 #扩胞情况
MP = 16 16 8 #Mesh网格密度
TMAX = 2000
TMIN = 0
TSTEP = 10 #温度上下限默认为100,0，step默认值为10

XYZ_PROJECTION = .TRUE.
打开该tag，将对特征向量按实空间进行投影，得到不同原子在a,b,c三个轴向上的位移投影，数据存储为frequency atom1_x atom1_y atom1_z atom2_x atom2_y atom2_z…形式。

PROJECTION_DIRECTION = 1 1 0
打开该tag，将原子位移在指定实空间轴向上进行投影。
phonopy -c POSCAR thermal_displacements.conf -p -s
执行后，输出结果将保存在thermal_displacements.yaml中。注意，计算的原子位移为均方位移。

此外的两个tag：

TDISPMAT = .TRUE.
打开该tag时，将以3*3矩阵形式储存原子在六个方向上的均方位移xx, yy, zz, yz, xz, xy，结果写入thermal_displacement_matrices.yaml中。

TDISPMAT_CIF = 1273.0
该tag用于指定计算在某温度（1273K）下，原子的均方位移矩阵，结果输出在tdispmat.cif中。

Phonopy其他功能

  • phonopy --symmetry -c POSCAR-unit

BPOSCAR和PPOSCAR生成，分别是标准胞和原胞，vasp格式。

  • phonopy --rd 10 --dim 2 2 2 --pa auto --amplitude 0.03 -c POSCAR-unitcell

–rd 10 代表生成10个随机位移结构，–dim 2 2 2扩胞倍数，–pa auto沿三个坐标轴投影方向，–amplitude 0.03位移幅度，-c指定对象。

输出”phonopy_disp.yaml” and supercells

  • phonopy-load phonopy_disp_orig.yaml --rd 1000 --temperature 300

generate 1000 supercells with displacements at 300K.

  • New in v2.25 Symfc (symfc/symfc) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use symfc, its python module has to be installed via conda-forge or pip.

New in v2.3 ALM (ttadano/ALM) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use ALM, its python module has to be installed via conda-forge or building it. The installation instruction is found here.

When ALM is used, please cite the paper: T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018).

FC_CALCULATOR = ALM