我的学习记录
我的学习记录利用phono3py在尝试讨论FGT声子与声子相互作用的练习中建立的
准备：
pw.FGT.in的输入文件：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
50
51
52
53
54
55
56
57
58&CONTROL
calculation='scf',
restart_mode='from_scratch',
wf_collect=.false.,
prefix='FGT',
pseudo_dir='./pseudo',
forc_conv_thr=1.0d-3
tddft_is_on=.True.
td_outputS=.True.
edt=3.0000
dt=3.0000
nstep=6202
diagonSteps=1
/
&SYSTEM
ibrav = 0
nat = 6
ntyp = 3
ecutwfc=120
nbnd=120
occupations='smearing'
smearing='gauss'
degauss=0.01
input_dft='pz'
starting_magnetization(1)=0
starting_magnetization(2)=0
starting_magnetization(3)=1
lspinorb=.TRUE.
nosym=.TRUE.
noinv = .True.
noncolin=.TRUE.
/
&ELECTRONS
electron_maxstep=200
conv_thr=1.0d-6
/
&IONS
/
&CELL
/
CELL_PARAMETERS {angstrom}
3.8564786911 0.0000000000 0.0000000000
-1.9340042258 3.3364748978 0.0000000000
0.0000000000 0.0000000000 17.4066505432
ATOMIC_SPECIES
Te 127.60000 Te.upf
Ge 72.630000 Ge.upf
Fe 55.84500 Fe.upf
ATOMIC_POSITIONS {crystal}
Te 0.666678011 0.333322078 0.893757522
Te 0.666677892 0.333322078 0.606242478
Ge 0.333640367 0.666359782 0.750000179
Fe 0.999609590 0.000390340 0.818844140
Fe 0.999609590 0.000390339 0.681155801
Fe 0.667117894 0.332882077 0.749999881
K_POINTS {automatic}
11 11 1 0 0 0
head.in:
首先基于pw.FGT.in产生具有位移的不同的超胞文件以及phono3py_disp.yaml：
1 | phono3py --qe -d --dim="2 2 1" -c pw.FGT.in |
这一个步骤中会生成多个超胞*.in文件
注意：这一步骤中，只能用以下格式ATOMIC_POSITIONS {crystal}且CELL_PARAMETERS {angstrom}：
1 | CELL_PARAMETERS {angstrom} |
其次，对于超胞计算力
利用脚本在disp-xxx中计算关于超胞supercell-xxx.in并在每一次计算后在当前目录下输出supercell-xxx.out1
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
50
51#!/bin/bash
#
#SBATCH --job-name=ph-ph
#SBATCH --output=FGT.out
#SBATCH -N 5
#SBATCH --ntasks-per-node=52
##SBATCH --time=3-00:00:00
#SBATCH -p long
ulimit -s unlimited
ulimit -c unlimited
#module load pmix/2.2.2
module load parallel_studio/2020.2.254
module load intelmpi/2020.2.254
#EXEC=/home/users/nawu/qe/TDPW6.6/pw.x
#EXEC=/home/users/nawu/spin/tdpw.x
#spin-direction
#EXEC=/home/users/nawu/TDPW/new/src/tdpw.x
EXEC=/home/users/nawu/new/tdpw.x
#EXEC=/home/users/nawu/software/TDPW6.6/pw.x
#EXEC=/home/users/nawu/qe/q-e-qe-6.6/bin/epw.x
mkdir supercell_out
for i in {00001..02760}
do
cp -r pseudo supercell_out
cat head.in supercell-$i.in >| new-$i.in
mv new-$i.in supercell_out
cd supercell_out
mkdir disp-$i
mv pseudo disp-$i/
mv new-$i.in disp-$i/
cd disp-$i
srun --mpi=pmi2 $EXEC -npools 5 -i new-$i.in | tee new-$i.out
cd
cd FGT/new/ph-ph
done
#srun --mpi=pmi2 $EXEC -npool 5 -i pw.FGT.in | tee pw.FGT.out
#srun --mpi=pmi2 $EXEC -i pw.MoS2.nscf-SOC.in | tee pw.MoS2.nscf-SOC.out
#srun --mpi=pmi2 $EXEC -i pw.MoS2.bands-SOC.in | tee pw.MoS2.bands-SOC.out
#srun --mpi=pmi2 $EXEC -i MoS2.band.in | tee MoS2.band.out
#mpirun $EXEC -i input.in | tee result
exit
然后，获得力
通过--cf3选项获得在QE计算中每个原子的受力情况，进一步获得力 FORCES_FC31
phono3py --cf3 disp-{00001..02760}/new-{00001..02760}.out
对于phono3py_disp.yaml是用来创建 FORCES_FC3，因此他必须存在于当前目录下
接着，计算二阶和三阶力常数
也就是创建 fc3.hdf5 和 fc2.hdf5文件
1
phono3py --sym-fc
这里--sym-fc是对于fc3和fc2对称化
最后，如果有需要，计算晶格热导（LTC— lattice thermal conductivity)、
1 | phono3py --mesh="11 11 11" --fc3 --fc2 --br |
■Theory
DFPT—-参见DFPT学习的笔记
■ WorkFlow
Phono3py计算的是声子-声子相互作用以及相关性质。
下图阐释了晶格热导（lattice thermal conductivity—-LTC）
其他的性质比如谱函数可以通过以下workflow来实现
LTC计算通过下面三个步骤来实现：
(1)计算超胞力常数集（force-sets）
(2)计算力常数（FC）
(3)计算晶格热导（lattice thermal conductivity—-LTC）
用户在每一步调用phono3py时，至少将超胞（supercell matrix:supercell+displacements)或单胞(unit cell)作为输入
第一步:
由产生的超胞和原子位置的集合（位移集合称之为”位移数据集”）
具有位移的超胞这里生成
这些超胞模型的力可以通过第一性原理计算得到一个”力集
第二步:
由第一步得到的位移集和力集来计算二阶（fc2）和三阶(fc3)力常数
第三步：
由第二步得到的力常数来计算晶格热导。当输入的不是一个初基晶胞，需要被给定初基晶胞。
当提供了非解析项修正的参数时，将会包含长程的dipole-dipole相互作用
■ 设置参数
phonopy的设置tag和命令可以同时使用，但如果同时使用，优先使用命令
例如：下面是一个配置案例
储存在setting.conf文件中1
2
3
4
5
6
7
8
9DIM = 2 2 2
DIM_FC2 = 4 4 4
PRIMITIVE_AXES = 0 1/2 1/2 1/2 0 1/2 1/2 1/2 0
MESH = 11 11 11
BTERTA = .TRUE.
NAC = .TRUE.
READ_FC2 = .TRUE.
READ_FC3 = .TRUE.
CELL_FILENAME = POSCAR-unitcell
这里的设置参数不区分大小写。
可以通过以下命令来运行：1
2
3phono3py setting.conf [comannd options]
或者
phono3py [comannd options] -- setting.conf
输入参数
⚪输入原胞文件名字
·-c (CELL_FILENAME)
1 | phono3py -c POSCAR-unitcell ... (many options) |
⚪计算接口
·—qe (CALCULATOR = QE)
涉及到qe
·—crystal (CALCULATOR = CRYSTAL)
涉及到crystal
·—turbomole (CALCULATOR = TURBOMOLE)
涉及到turbomole
⚪创建默认输入参数的用法
·—cf3 (command option only)
这用于从phono3py_disp.yaml和包含超单元中力的力计算中输出创建FORCES_FC3
phono3py_disp.yaml 位于当前目录中。必须指定计算器接口，但 VASP（默认）情况除外。1
2phono3py --cf3 disp-{00001..00755}/vasprun.xml
phono3py --cf3 supercell_out/disp-{00001..00111}/Si.out
·—cf3-file (command option only)
1 | phono3py --cf3-file file_list.dat |
file_list.dat包含可以从当前目录中识别的文件名，并且应如下所示：1
2
3
4
5disp-00001/vasprun.xml
disp-00002/vasprun.xml
disp-00003/vasprun.xml
disp-00004/vasprun.xml
...
文件名的顺序很重要。此选项与--cutoff-pair 一起使用
·—cf2 (command option only)
这用于创建类似于 —cf3 选项FORCES_FC2。
phono3py_disp.yaml 必须位于当前目录中。这是可选的。
必须指定计算接口，但 VASP（默认）情况除外。
FORCES_FC2必须使用 —dim-fc2 选项运行
1 | phono3py --cf2 disp_fc2-{00001..00002}/vasprun.xml |
·—cfz (command option only)
这与 —cf3 和 —cf2 组合,用于创建FORCES_FC3和FORCES_FC2的残余力。
必须指定计算接口，但 VASP（默认）情况除外。
在下面的示例中，假设 disp3-00000/vasprun.xml 和 disp2-00000/vasprun.xml包含完美超胞的力。在理想情况下，这些力为零，但通常不是。在这里，这被称为”残余力（residual forces）”。
有时，力常数的质量以这种方式得到改善
1
2phono3py --cf3 disp3-{00001..01254}/vasprun.xml --cfz disp3-00000/vasprun.xml
phono3py --cf2 disp2-{00001..00006}/vasprun.xml --cfz disp2-00000/vasprun.xml
·—fs2f2 or —force-sets-to-forces-fc2 (command option only)
FORCES_FC2是从phonopy的FORCE_SETS文件创建的。
phono3py_disp.yaml 所需的 yaml 行显示为文本。1
phono3py --fs2f2
与 –dim-fc2 结合使用时，phonopy 的FORCE_SETS是从 FORCES_FC2 和 phono3py_disp.yaml 创建的，而不是 FORCES_FC3 和 phono3py_disp.yaml。1
phono3py --cfs --dim-fc2="x x x"
·—cfs or —create-force-sets (command option only)
Phonopy的FORCE_SETS是从FORCES_FC3和phono3py_disp.yaml创建的。1
phono3py --cfs
⚪超胞和原胞
·—dim (DIM)
指定超胞维度，如果当前目录下包含phono3py_disp.yaml，不必要指定
·—dim-fc2 (DIM_FC2)
指定二阶力常数（对于简谐声子）的超胞维度。这是可选的。
当前目录中存在正确的phono3py_disp.yaml 时，无需指定。
对于二阶力常数，可以使用此选项指定比三阶力常数更大且不同的超胞大小。
通常，一对原子之间的相互作用在实际空间中比三个原子之间的相互作用具有更长的范围。
因此，为了减少计算需求，仅为二阶力常数选择更大的超单元尺寸可能是一个好主意
通过利用-d选项，结构文件和phono3py_disp.yaml被创建
这些被用来计算较大超胞中的二阶力常数，也可以计算三阶力常数1
phono3py -d --dim="2 2 2" --dim-fc2="4 4 4" -c POSCAR-unitcell
在力计算后，—cf2选项被用来创建FORCES_FC21
phono3py --cf2 disp-{001,002}/vasprun.xml
为了计算较大超胞的二阶力常数，必须要准备：FORCES_FC2， phono3py_disp.yaml
当运行 phono3py计算较大的二阶力常数时，必须指定 --dim-fc2选项。
（由于运行 phono3py 而创建的 fc2.hdf5 包含具有较大超胞的二阶力常数。）
文件名与在没有 --dim-fc2 选项的常规 phono3py 运行中创建的文件名相同。1
phono3py --dim="2 2 2" --dim_fc2="4 4 4" -c POSCAR-unitcell
·—pa, —primitive-axes (PRIMITIVE_AXES)
从非初基原胞变换到初基晶胞
⚪创建位移
·-d (CREATE_DISPLACEMENTS = .TRUE.)
创建具有位移的超胞。
使用 --amplitude选项，可以控制原子位移距离。
使用此选项，将创建具有位移的超胞的文件和phono3py_disp.yaml 文件。
如果输入单胞结构不是初基原胞，则应指定 --pa，例如，如果输入单元单元具有 F 中心，则应指定 --pa="F"。
·—amplitude (DISPLACEMENT_DISTANCE)
指定原子位移的距离
⚪创建力常数
·—cfc or —compact-fc (COMPACT_FC = .TRUE.)
当通过 FORCES_FC3 and/or FORCES_FC2,使用较小尺寸的力常数会被产生
fc2的数据形式：(num_patom, num_satom)
fc3的数据形式：(num_patom, num_satom, num_satom)
这里num_patom是指原胞中的原子数,num_satom是指超胞中的原子数
在全尺寸力常数的情况下，num_patom被num_satom所取代。1
phono3py --dim="2 2 2" --cfc --pa="F" -c POSCAR-unitcell
·—sym-fc (FC_SYMMETRY = .TRUE.)
二阶力常数和三阶力常数对称性可以通过-sym-fc2 (SYMMETRIZE_FC2 = .TRUE.)
和 —sym-fc3r (SYMMETRIZE_FC3 = .TRUE.)分别应用
·—cutoff-fc3 or —cutoff-fc3-distance (CUTOFF_FC3_DISTANCE)
此选项不用于减少具有位移的超胞的数量，
但此选项用于在给定三阶力常数的元素中设置零。
零元素的选择条件是每个原子三元组中原子的任何成对距离大于指定的截止距离。
如果要减少超胞的数量，第一个选择是减小超胞的大小，第二个选择是使用—cutoff-pair选项
·—cutoff-pair or —cutoff-pair-distance (CUTOFF_PAIR_DISTANCE)
⚪倒空间取样
·—mesh (MESH or MESH_NUMBERS)
在倒空间里通过指定数目生成mesh取样网格，这样的取样网格总是以Gamma点为中心沿着倒格矢
·—gp (GRID_POINTS)
网格由唯一索引指定，网格点的索引由空格或逗号 （，） 分隔的数字指定。
网格点与其索引之间的映射表是通过运行 —loglevel=2 选项获得的。1
phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --write-gamma --gp="0 1 2 3 4 5"
这里--gp="0 1 2 3 4 5"可以写成--gp=0,1,2,3,4,5
·—ga (GRID_ADDRESSES)
这用于指定网格点，如 —gp 选项，但其地址由整数表示。
·—bi (BAND_INDICES)
指定能带指数，以空格分隔的索引处的计算值取平均值，并单独计算以逗号分隔的值。
输出文件名将是，例如，gammas-mxxx-gxx（-sx）-bx.dat其中 bxbx…显示了波段指数，其中值在这些波段上进行计算、求和和和平均1
phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --gp="34" --bi="4 5, 6"
·—wgp (command option only)
·—stp (command option only)
⚪布里渊区积分
·—thm (TETRAHEDRON = .TRUE.)
四面体法用于计算自能量的虚部。这是默认选项。
因此，除非预期一次性执行中通过tetrahedral方法和smearing方法获得的结果，否则不必指定此值
·—sigma (SIGMA)
对于smearing中计算自能虚部的σ的数值
多个σ值也由空格分隔的数值指定。
·—sigma-cutoff (SIGMA_CUTOFF_WIDTH)
指定数值精确度
·—full-pp (FULL_PP = .TRUE.)
⚪求解声子Boltzman方程
·—br (BTERTA = .TRUE.)
·—lbte (LBTE = .TRUE.)
⚪散射
·—isotope (ISOTOPE =.TRUE.)
·—mass-variances or —mv (MASS_VARIANCES)
·—boundary-mfp, —bmfp (BOUNDARY_MFP)
·—ave-pp (USE_AVE_PP = .TRUE.)
·##### ·—const-ave-pp (CONST_AVE_PP = .TRUE.)
·—nu (N_U = .TRUE.)
·—scattering-event-class (SCATTERING_EVENT_CLASS)
Temperature
·—ts (TEMPERATURES): Temperatures
·—tmax, —tmin, —tstep (TMAX, TMIN, TSTEP)
Non-analytical term correction
·—nac (NAC = .TRUE.)
·—q-direction (Q_DIRECTION)
Imaginary and real parts of self energy
·—ise (IMAG_SELF_ENERGY = .TRUE.)
·—rse (REAL_SELF_ENERGY = .TRUE.)
Spectral function
·—spf (SPECTRAL_FUNCTION = .TRUE.)
Joint density of states (JDOS) and weighted-JDOS
·—jdos (JOINT_DOS = .TRUE.)
Sampling frequency for distribution functions
·—num-freq-points, —freq-pitch (NUM_FREQUENCY_POINTS)
Mode-Gruneisen parameter from 3rd order force constants
·—gruneisen (GRUNEISEN = .TRUE.)
File I/O
·—fc2 (READ_FC2 = .TRUE.)
从fc2.hdf5中读取二阶力常数
·—fc3 (READ_FC3 = .TRUE.)
从fc3.hdf5中读取三阶力常数
·—write-gamma (WRITE_GAMMA = .TRUE.)
简谐声子频率下自能量的虚构部分以hdf5格式写入文件。
结果被写入 kappa-mxxx-gx（-sx-sdx）.hdf5 或 kappa-mxxx-gx-bx（-sx-sdx）.hdf5 与 —bi 选项。使用 —sigma 和 —sigma-cutoff 选项，-sx 和 —sdx 分别插入到 .hdf5 的前面。
·—read-gamma (READ_GAMMA = .TRUE.)
简谐声子频率下的自能量的虚部以hdf5格式从kappa文件中读取。
最初搜索 kappa-mxxx（-sx-sdx）.hdf5 的常用结果文件。
除非找到它，否则它会尝试读取每个网格点的 kappa 文件 kappa-mxxx-gx（-sx-sdx）.hdf5。然后，类似地，找不到 kappa-mxxx-gx（-sx-sdx）.hdf5，搜索 kappa-mxxx-gx-bx（-sx-sdx）.hdf5 文件的波段索引。
·—write-gamma-detail (WRITE_GAMMA_DETAIL = .TRUE.)
每个q点三元组分对自能量的虚部的贡献都写入gamma_detail-mxxx-gx（-sx-sdx）.hdf5文件中。
·—write-phonon (WRITE_PHONON = .TRUE.)
声子频率、特征向量和网格点地址存储在声子-mxxx.hdf5 文件中。可能需要 –pa 和 –nac，具体取决于计算设置。1
phono3py --fc2 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --write-phonon
·—read-phonon (READ_PHONON = .TRUE.)
从 phonon-mxxx.hdf5 文件中读取声子频率、特征向量和网格点地址，并使用这些声子值继续计算。1
phono3py --fc2 --fc3 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --read-phoonon --br
·—write-pp (WRITE_PP = .TRUE.) and —read-pp (READ_PP = .TRUE.)
声子和声子(ph-ph)之间的相互作用强度从pp-mxxx-gx.hdf5中写入和读取
这仅适用于晶格热导率的计算，即只能与 --br 或--lbte 一起使用
使用和不使用指定--full-pp选项时，存储的数据是不同的
在前一种情况下，所考虑的声子三元组之间的所有ph-ph相互作用强度都以简单的方式存储，
但在后一种情况下，只有计算碰撞的必要元素以复杂的方式存储。1
2phono3py --fc2 --fc3 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --write-pp --br --gp=1
phono3py --fc2 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --read-pp --br --gp=1
·—hdf5-compression (command option only)
默认情况下，大多数 phono3py HDF5 输出文件都使用 gzip 压缩过滤器进行压缩。
为避免压缩，必须设置 —hdf5-compression=None。
·-o (command option only)
修饰默认输出的文件名，例如,-o iso 文件名kappa-m191919.hdf5被更改为kappa-m191919.iso.hdf5
这个规则适用于以下文件：
fc3.hdf5
fc2.hdf5
kappa-xxx.hdf5
phonon-xxx.hdf5
pp-xxx.hdf5
gamma_detail-xxx.hdf5 (write only)
·-i (command option only)
修饰默认的输入文件，例如指定-i iso —fc3，文件名是fc3.iso.hdf5将会代替fc3.hdf5
这个规则适用于以下文件：
fc3.hdf5
fc2.hdf5
kappa-xxx.hdf5
phonon-xxx.hdf5
pp-xxx.hdf5
·—io (command option only)
同时修饰输入和输出文件
■ 输出文件
（1）中间的文本文件
以下文件不与phonopy兼容，但是phonopy创立的FORCE_SETS文件可以通过phono3py的命令创建
·phono3py_disp.yaml
可通过-d选项创建
·FORCES_FC3
可通过—cf3选项创建，有两种类型,第一种类型:
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 # File: 1
# 1 0.0300000000000000 0.0000000000000000 0.0000000000000000
-0.6458483000 0.0223064300 -0.0143299700
0.0793497000 0.0088413200 -0.0052766800
0.0768176500 -0.0095501600 0.0057262300
-0.0016552800 -0.0366684600 -0.0059480700
-0.0023432300 0.0373490000 0.0059468600
0.0143901800 0.0000959800 -0.0001100900
-0.0019487200 -0.0553591300 -0.0113649500
0.0143732700 -0.0000614400 0.0000502600
-0.0020311400 0.0554678300 0.0115355100
...
# File: 1254
# 37 0.0000000000000000 0.0000000000000000 -0.0300000000000000
# 68 0.0000000000000000 0.0000000000000000 -0.0300000000000000
-0.0008300600 -0.0004792400 0.0515596200
-0.0133197900 -0.0078480800 0.0298334900
0.0141518600 -0.0105405200 0.0106313000
0.0153762500 -0.0072671600 0.0112864200
-0.0134565300 -0.0076112400 0.0298334900
-0.0019180000 -0.0011073600 0.0272454300
0.0013945800 0.0169498000 0.0112864200
0.0006578200 0.0003797900 0.0085617600
-0.0020524300 0.0175261300 0.0106313000
0.0019515200 0.0011267100 -0.2083651200
0.0148675800 -0.0516285500 -0.0924200600
-0.0168043800 0.0074232400 -0.0122506300
-0.0128831200 0.0114004400 -0.0110906700
...
这个文件包含超胞的力，以#开始的部分可以忽略，每一行给出原子所受到的以笛卡尔坐标表示的力
力都储存在一个numpy数组中（forces）以(num_supercells, num_atoms_in_supercell, 3)格式
第二类格式同phonopy的第Ⅱ类格式
·FORCES_FC2
可通过—cf2选项创建
（2）HDF5文件
·kappa-*.hdf5
·fc3.hdf5
（以eV/Angstrom3为单位储存起来的实空间的三阶力常数矩阵）
在phono3py中，它以 dtype='double' 和 order='C'的numpy矩阵以下列形式储存起来:
1
(num_atom, num_atom, num_atom, 3, 3, 3)反应的是`Φ(α,β,γ)(lk,l'k',l''k'')`
前三个num_atom是代表着超胞中对应lk,l'k',l''k''的原子索引，后三个元素代表着对应α,β,γ的笛卡尔坐标
如果想输入一个超胞结构和fc3，需要匹配单胞和超胞之间的原子索引
（这看起来很容易的由phono3py将超胞当作单胞来处理，如POSCAR,unitcell.in,且指定初基晶胞可以用-pa的选项）
例如：考虑一个超胞扩胞至单胞没有中心的222大小的，那么--pa参数应该设置为 1/2 0 0 0 1/2 0 0 0 1/2
如果单胞是一个惯用晶胞且由中心，如面心，那么--pa参数应该设置为 0 1/4 1/4 1/4 0 1/4 1/4 1/4 0
·fc2.hdf5
（以eV/Angstrom2为单位储存起来的实空间的二阶力常数矩阵）
在phono3py中，它以 dtype='double' 和 order='C'的numpy矩阵以下列形式储存起来:
1
(num_atom, num_atom, 3, 3)反应的是`Φ(α,β)(lk,l'k')`
·gamma-*.hdf5
以THz为单位的谐振声子频率下自能的虚部（ Γλ(ωλ)= 半线宽）
·gamma_detail-*.hdf5
以THz为单位的谐振声子频率下对于自能的虚部Q点的三重贡献
（3）简单的文本文件
·gammas-*.dat
以THz为单位对应频率为Γλ(ωλ)自能的虚部
·jdos-*.dat
以THz为单位的联合DOS（JDOS）
·linewidth-*.dat
Notes：
计算结果的输出文件大多储存在HDF5格式的文件中
■ 如何读取储存在hdf5文件中的结果？
@ https://docs.h5py.org/en/latest/high/dataset.html#reading-writing-data%3E
How to use HDF5 python library
1 | "MgO的热导率计算结果以及在300K上看观察热导率" |
Details of kappa-*.hdf5 file
文件名，例如 kappa-m323220.hdf5，由某些特定选项确定。
mxxx，显示采样网格的数量。
sxxx 和 gxxx 可选择出现。
sxxx 给出了关于在布里渊积分中计算声子寿命的smearing法中的smearing宽度，gxxx 表示网格数目。
使用命令选项 -o，可以稍微修改文件名。例如，-o nac给出了kappa-m323220.nac.hdf5来记住选项—nac被使用。
目前，kappa-*.hdf5 文件（不适用于特定的网格点）包含如下所示的属性。
mesh
沿着倒格矢的网格取样点
frequency
声子的频率以THz为单位，矩阵形状（不可约表示的q点，声子谱）
gamma
自能的虚部，单位是THz，
矩阵形状是：
·针对所有网格点的（温度，不可约表示的q点，声子谱）
·针对特定网格点的（温度，声子谱）
声子的寿命可以通过gamma来估计
2pi必须与 hdf5 文件中的 gamma 值相乘，才能将频率单位转换为角频率。
gamma_isotope
group_velocity
声子群速度，物理单位是THz
阵列形状为（不可约表示的q点，声子谱，3=笛卡尔坐标）
heat_capacity
比热容，单位：eV/K
阵列形状为（温度，不可约表示的q点，声子谱）
kappa
热导张量，单位：W/m-K
矩阵形状为（温度，6=（xx,yy,zz,yz,xz,xy））
mode-kappa
在k*处的热导张量。
gv_by_gv
q-point
在约化的坐标中不可约表示的q点
矩阵形状：（不可约表示的q点，3=倒空间的约化坐标）
temperature
这里计算热导率时的温度,单位为K
weight
对应不可约表示的q点上的权重，权重求和等于网格点数目
ave_pp
平均的声子与声子间相互作用，单位eV2
kappa_unit_conversion
How to know grid point index number corresponding to grid address
■辅助工具
phono3py-kaccum（晶格热导和相关性质）
phono3py-kdeplot（声子寿命的分布）
■如何使用
以Si-example为例
为了避免二次计算fc3和fc2，fc3.hdf5和fc2.hdf5在单节点上进行计算1
phono3py --dim="2 2 2" --sym-fc -c POSCAR-unitcell
不可约网格点的索引需要由 –wgp 选项找到指定 —ga 选项1
phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --wgp
且储存在ir_grid_points.yaml1
2egrep '^- grid_point:' ir_grid_points.yaml|awk '{printf("%d,",$3)}'
0,1,2,3,4,5,6,7,8,9,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,60,61,62,63,64,65,66,67,68,69,70,71,72,73,80,81,82,83,84,85,86,87,88,89,90,91,100,101,102,103,104,105,106,107,108,109,120,121,122,123,124,125,126,127,140,141,142,143,144,145,160,161,162,163,180,181,402,403,404,405,406,407,408,409,422,423,424,425,426,427,428,429,430,431,432,433,434,435,442,443,444,445,446,447,448,449,450,451,452,453,462,463,464,465,466,467,468,469,470,471,482,483,484,485,486,487,488,489,502,503,504,505,506,507,522,523,524,525,542,543,804,805,806,807,808,809,824,825,826,827,828,829,830,831,832,833,844,845,846,847,848,849,850,851,864,865,866,867,868,869,884,885,886,887,904,905,1206,1207,1208,1209,1226,1227,1228,1229,1230,1231,1246,1247,1248,1249,1266,1267,1608,1609,1628,1629,
上面显示的所有网格点的计算数据作为指数是获得晶格热导率所必需的。
为了将计算需求分配到计算机节点，选择并执行一组网格点索引，如下所示：1
phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --gp="0,1,2,3,4,5,6,7,8,9,20,21,22,23,24,25" --write-gamma
然后生成许多kappa-m191919-gx.hdf5文件。这些文件名不应该被更改，因为在通过phono3py读取数据时，这些文件名应该是这样，尽管有一点自由来排列这些文件名，为此请参阅-o和-i选项。完成所有不可约网格点指数的计算后，RTA热导率将在短时间内通过存储的数据进行另一次运行计算：1
phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --read-gamma
■使用cutoff pair-distance计算力常数
使用 --cutoff pair-distance，可减少具有要计算的位移的超胞的数量。
当然，这牺牲了三阶力常数（fc3）的准确性。
在phono3py中默认计算超胞中所有的三阶力常数，尽管这很依赖于超胞中原子的数目和晶体的对称性，原子对的状态很大已经超过了我们的计算资源
有时，我们期盼fc3在三个原子之间的相互作用范围比选择的超胞的大小小一点，因此我们需要略去一部分计算的元素，这可以利用这个选项来实现
一个超胞的fc3元素可以通过三个原子位移来指定，三个中的两个Δr1和Δr2是有限位移的（和），但其中一个包含在给出力的计算中
截断半径dcut可以定义为两个原子之间1和2距离的上边界，由此，原子对之间的距离集合{（r1,r2）||r1-r2|<dcut}来确定fc3的计算
如果三个原子之间的距离都大于dcut.那么fc3元不会包含在内，而是会简单的设置为0
用法
创建具有位移的超胞
--cutoff-pair选项是用来创建具有位移的超胞的，因此这个选项通常结合-d来使用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
26hono3py --cutoff-pair=5 -d --dim="2 2 2" -c POSCAR-unitcell
...
Unit cell was read from "POSCAR-unitcell".
Displacement distance: 0.03
Number of displacements: 111
Cutoff distance for displacements: 5.0
Number of displacement supercell files created: 51
Displacement dataset was written in "phono3py_disp.yaml".
...
% ls POSCAR-0*
POSCAR-00001 POSCAR-00032 POSCAR-00043 POSCAR-00080 POSCAR-00097
POSCAR-00002 POSCAR-00033 POSCAR-00070 POSCAR-00081 POSCAR-00098
POSCAR-00003 POSCAR-00034 POSCAR-00071 POSCAR-00082 POSCAR-00099
POSCAR-00016 POSCAR-00035 POSCAR-00072 POSCAR-00083 POSCAR-00100
POSCAR-00017 POSCAR-00036 POSCAR-00073 POSCAR-00084 POSCAR-00101
POSCAR-00018 POSCAR-00037 POSCAR-00074 POSCAR-00085 POSCAR-00102
POSCAR-00019 POSCAR-00038 POSCAR-00075 POSCAR-00086 POSCAR-00103
POSCAR-00024 POSCAR-00039 POSCAR-00076 POSCAR-00087
POSCAR-00025 POSCAR-00040 POSCAR-00077 POSCAR-00088
POSCAR-00026 POSCAR-00041 POSCAR-00078 POSCAR-00089
POSCAR-00027 POSCAR-00042 POSCAR-00079 POSCAR-00096
% ls POSCAR-0*|wc -l
51
总结:
Number of displacements: 111表示的是：
如果没有--cutoff-pair选项，具有位移的超胞的数目
displacement supercell files created: 51表示的是：
通过 —cutoff-pair选项给出具有位移的超胞的收缩数量，这里POSCAR数目是51
这一步骤中phono3py_disp.yaml被创建，这包含有关此收缩的信息，并在其他计算步骤中使用，因此必须仔细保存此文件。
超胞文件
如果要置换的一对原子之间的距离大于指定的截止对距离，则不会生成POSCAR-xxxxx（在其他计算软件中，文件名的前缀不同）
特殊的phono3py_disp.yaml
将 --cutoff-pair选项与 -d 选项一起使用，将创建一个特殊的 phono3py_disp.yaml。
这包含有关位移原子对之间的距离以及是否要计算这些原子对的信息。
此特殊phono3py_disp.yaml对于创建 fc3 是必需的，因此请注意不要通过运行不带 --cutoff-pair或具有不同值的不同--cutoff-pair的选项-d来覆盖它。
创建FORCES_FC3
要创建FORCES_FC3，只有使用 --cutoff-pair选项创建的超胞的输出文件才会作为参数传递给 phono3py。
特殊的phono3py_disp.yaml 文件必须位于当前目录中。
Si的示例：
在这里，假设力是使用disp-xxxxx目录中的VASP计算的。
运行力计算后，每个目录中都应该有包含力的输出文件（对于 VASP vasprun.xml）。1
2
3
4
5
6phono3py --cf3 disp-{00001,00002,00003,00016,00017,00018,00019,00024,00025,00026,00027,00032,00033,00034,00035,00036,00037,00038,00039,00040,00041,00042,00043,00070,00071,00072,00073,00074,00075,00076,00077,00078,00079,00080,00081,00082,00083,00084,00085,00086,00087,00088,00089,00096,00097,00098,00099,00100,00101,00102,00103}/vasprun.xml
...
Displacement dataset is read from phono3py_disp.yaml.
counter (file index): 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 50 51
FORCES_FC3 has been create
当力文件过大使用-cf3-file选项是推荐的。1
2for i in `ls POSCAR-0*|sed s/POSCAR-//`;do echo disp-$i/vasprun.xml;done > file_list.dat
% phono3py --cf3-file file_list.dat
计算声子-声子相互作用
为了创建 fc3，--cutoff-pair选项不是必需的，但需要特殊的 phono3py_disp.yaml
一旦创建了fc3.hdf5和fc2.hdf5,--cutoff-pair选项和特殊的phono3py_disp.yaml不再会被需要1
2
3
4
5phono3py --mesh="11 11 11" --fc3 --fc2 --br
...
300.0 108.051 108.051 108.051 0.000 0.000 -0.000
...
■ Phono3py API
·Crystal structure
phono3py中的晶体结构通常是PhonopyAtoms类实例。
当我们想从以力计算器格式写入的文件中获取PhonopyAtoms类实例时，
可以使用phonopy中的read_crystal_structure函数，例如，在AlN-LDA示例中，1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22In [1]: from phonopy.interface.calculator import read_crystal_structure
In [2]: unitcell, _ = read_crystal_structure("POSCAR-unitcell", interface_mode='vasp')
In [3]: print(unitcell)
lattice:
- [ 3.110999999491908, 0.000000000000000, 0.000000000000000 ] # a
- [ -1.555499999745954, 2.694205030733368, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 4.978000000000000 ] # c
points:
- symbol: Al # 1
coordinates: [ 0.333333333333333, 0.666666666666667, 0.000948820000000 ]
mass: 26.981539
- symbol: Al # 2
coordinates: [ 0.666666666666667, 0.333333333333333, 0.500948820000000 ]
mass: 26.981539
- symbol: N # 3
coordinates: [ 0.333333333333333, 0.666666666666667, 0.619051180000000 ]
mass: 14.006700
- symbol: N # 4
coordinates: [ 0.666666666666667, 0.333333333333333, 0.119051180000000 ]
mass: 14.006700
或者，它是直接从PhonopyAtoms类创建的。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17In [1]: from phonopy.structure.atoms import PhonopyAtoms
In [2]: a = 3.11
In [3]: a = 3.111
In [4]: c = 4.978
In [5]: lattice = [[a, 0, 0], [-a / 2, a * np.sqrt(3) / 2, 0], [0, 0, c]]
In [6]: x = 1. / 3
In [7]: points = [[x, 2 * x, 0], [x * 2, x, 0.5], [x, 2 * x, 0.1181], [2 * x, x, 0.6181]]
In [8]: symbols = ['Al', 'Al', 'N', 'N']
In [9]: unitcell = PhonopyAtoms(cell=lattice, scaled_positions=points, symbols=symbols)
·Phono3py class
为了从python上使用phono3py，这里有一个phono3py API1
from phono3py import Phono3py
正如Workflow里写到的，phono3py的workflow可粗略的分为三步。
Phono3py类可以被利用到每一步,实例中的Ponon3py的最小输入参数是1
unitcell(第一个变量),supercell_matrix(超胞变量),primitive_matrix
他们这样使用：1
ph3 = Phono3py(unitcell, supercell_matrix=[3, 3, 2], primitive_matrix='auto')
unitcell是Phonopy class的例子，supercell_matrix 和primitive_matrix是通过unitcell来生成超胞和初基晶胞
这一步类似于Phononpy class
·Displacement dataset generation
Workflow中的(1)在超胞中产生位移，具有位移的超胞用作计算力的输入晶体结构模型
通过计算力，从phono3py中获得具有位移的超胞的力集，即phono3py仅提供晶体结构
这里,位移的集合称为displacement dataset,超胞的力的集合称为force sets,位移集和力集对作为简单的数据集，用于计算二阶和三阶力常数。
在安装完具有unitcell的phono3py后，如下位移产生：1
2
3
4
5
6
7
8
9In [4]: ph3 = Phono3py(unitcell, supercell_matrix=[3, 3, 2], primitive_matrix='auto')
In [5]: ph3.generate_displacements()
In [6]: len(ph3.supercells_with_displacements)
Out[6]: 1254
In [7]: type(ph3.supercells_with_displacements[0])
Out[7]: phonopy.structure.atoms.PhonopyAtoms
(由此，1254个具有位移的超胞被产生)
产生的位移集被用于力常数的计算中。
位移集和晶体结构信息被Phono3py.save()保存在一个文件中：1
In [8]: ph3.save("phono3py_disp.yaml")
·Supercell force calculation
生成的具有位移的超胞的力由一些第一性计算软件来计算
计算出的超胞的力将通过 Phono3py.forces属性存储在Phono3py类实例中，方法是设置形状为 （num_supercells， num_atoms_in_supercell， 3）的array_like变量。
在上面的示例中，数组形状为 （1254， 72， 3）。
如果计算的力集储存在FORCES_FC3文件中，forces的矩阵可由如下方式获得：1
2forces = np.loadtxt("FORCES_FC3").reshape(-1, num_atoms_in_supercell, 3)
assert len(forces) == num_supercells
·Force constants calculation
力常数的计算是通过位移集和力集对来实现的。phono3py.load从文件加载这些数据并在Phono3py类实例中设置它们的便捷功能。
但是，在仅预期为位移数据集的情况下，Phono3pyYaml 中的低级phono3py-yaml很有用。1
2
3
4
5
6
7In [1]: from phono3py.interface.phono3py_yaml import Phono3pyYaml
In [2]: ph3yml = Phono3pyYaml()
In [3]: ph3yml.read("phono3py_disp.yaml")
In [4]: disp_dataset = ph3yml.dataset
有了这个ph3yml，解释了如何计算力常数。在下文中，假定我们在当前目录中有FORCES_FC3。
位移数据集和力集设置为 Phono3py类实例，如下所示：1
2
3
4
5
6
7
8
9
10
11
12
13In [5]: unitcell = ph3yml.unitcell
In [6]: from phono3py import Phono3py
In [7]: ph3 = Phono3py(unitcell, supercell_matrix=ph3yml.supercell_matrix, primitive_matrix=ph3yml.primitive_matrix)
In [8]: import numpy as np
In [9]: forces = np.loadtxt("FORCES_FC3").reshape(-1, len(ph3.supercell), 3)
In [10]: ph3.dataset = disp_dataset
In [11]: ph3.forces = forces
由上，很容易的计算力常数：1
In [12]: ph3.produce_fc3()
·Non-analytical term correction parameters
用户从实验或计算中收集Born有效电荷和介电常数张量，并用作非解析项校正（NAC）的参数。这些参数存储为 python 字典：1
2
3
4
5
6
7
8'born': ndarray
Born effective charges
shape=(num_atoms_in_primitive_cell, 3, 3), dtype='double', order='C'
'factor': float
Unit conversion factor
'dielectric': ndarray
Dielectric constant tensor
shape=(3, 3), dtype='double', order='C'
请注意，在使用phono3py API的情况下，原始细胞中所有原子的Born有效电荷是必要的，而在BORN文件中，仅写入对称独立原子的Born有效电荷。
此 NAC 参数通过 Phono3py.nac_params属性设置为 Phono3py 类实例。可以从BORN 文件中读取NAC 参数1
2
3
4
5
6In [13]: from phonopy.file_IO import parse_BORN
In [14]: nac_params = parse_BORN(ph3.primitive, filename="BORN")
In [15]: ph3.nac_params = nan_params
其中，parse_BORN函数需要PhonopyAtoms类的初基原胞
·Regular grid for q-point sampling
声子通常取样在倒空间里以Gamma点为中心的一般网格。
这里有三种方式指定一般的网格：
1）Three integer values
2）One value
3）3x3 integer matrix
其中一个设置为mesh_numbers Phono3py类的属性。例如1
2
3ph3.mesh_numbers = [10, 10, 10]
ph3.mesh_numbers = 50
ph3.mesh_numbers = [[-10, 10, 10], [10, -10, 10], [10, 10, -10]]
·Three integer values
指定三个整数 (n1,n2,n3) 传统的常规网格由三个整数值定义。
q点由倒空间的基底与对应的整数线性组合给出：
·One value
指定一个距离值l
常规网格的三个整数值由以下计算定义：
此外，还支持使用广义正则网格。
通过在Phono3py 类实例化中指定 use_grg = True，将使用值l生成广义正则网格。
首先，搜索初基晶胞的惯用晶胞。
其次，将方程（2）应用于惯用晶胞的倒空间基底中。
然后，根据方程（1）相对于惯用晶胞的倒空间基底对q点进行采样。
·3x3 integer matrix (experimental)
他用于显式定义广义正则网格。广义正则网格被设计为由一个给定值自动生成，同时考虑对称性。但是，如果 3x3 整数矩阵（数组）正确遵循对称性，则接受此矩阵
·Phonon-phonon interaction calculation
在定义正则网格后计算三个声子相互作用强度。
当我们有上面使用的 phono3py_disp.yaml和FORCES_FC3（也可能是BORN）时，获取带有设置的Phono3py类实例很方便。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21In [1]: import phono3py
In [2]: ph3 = phono3py.load("phono3py_disp.yaml", log_level=1)
NAC params were read from "BORN".
Displacement dataset for fc3 was read from "phono3py_disp.yaml".
Sets of supercell forces were read from "FORCES_FC3".
Computing fc3[ 1, x, x ] using numpy.linalg.pinv with displacements:
[ 0.0300 0.0000 0.0000]
[ 0.0000 0.0000 0.0300]
[ 0.0000 0.0000 -0.0300]
Computing fc3[ 37, x, x ] using numpy.linalg.pinv with displacements:
[ 0.0300 0.0000 0.0000]
[ 0.0000 0.0000 0.0300]
[ 0.0000 0.0000 -0.0300]
Expanding fc3.
fc3 was symmetrized.
Max drift of fc3: 0.000000 (yxx) 0.000000 (xyx) 0.000000 (xxy)
Displacement dataset for fc2 was read from "phono3py_disp.yaml".
Sets of supercell forces were read from "FORCES_FC3".
fc2 was symmetrized.
Max drift of fc2: 0.000000 (yy) 0.000000 (yy)
声子之间的相互作用通过以下命令运行：1
2
3In [3]: ph3.mesh_numbers = 30
In [4]: ph3.init_phph_interaction()
三个声子之间的相互作用计算在 phono3py.phonon3.interaction.Interaction类中实现。
通过phono3py的设计选择，给出具有确定的q1下，遍取q2，三个声子之间(q1,q2,q3)的相互作用。
采用对称性是为了避免在固定q1处重复计算q2
一个Interaction储存固定q1处的数据。
三个声子之间的相互作用，需要较大的内存空间。
因此，在晶格导热系数计算中，它是按需计算的，然后在使用后放弃数据。
对于常规应用，三声子相互作用强度计算是计算要求最高的部分。详细的技术用于避免三声子相互作用强度的元素，如果它被指定这样做。
因此需要正确配置计算资源！
·Lattice thermal conductivity calculation
一旦ph3.init_phph_interaction（）计算了三个声子之间的相互作用，就可以计算晶格热导。
文档字符串中解释了许多参数，但即使使用默认参数，在弛豫时间近似模式下的晶格热导率计算也按如下方式执行：1
In [5]: ph3.run_thermal_conductivity()
·Use of different supercell dimensions for 2nd and 3rd order FCs
Phono3py支持二阶和三阶力常数（fc2和fc3）的不同超胞尺寸。
默认设置是使用相同的维度。
尽管 fc3 需要比 fc2 多得多的超级单元，但根据我们的经验，fc3 在直接空间中的相互作用范围更短。
因此，我们倾向于期望对 fc3 使用比 fc2 更小的超胞维度。
为了实现这一点，存在phonon_supercell_matrix参数来独立指定 fc2 超胞尺寸。1
2
3
4
5
6
7
8
9
10
11In [1]: from phonopy.interface.calculator import read_crystal_structure
In [2]: unitcell, _ = read_crystal_structure("POSCAR-unitcell", interface_mode='vasp')
In [3]: from phono3py import Phono3py
In [4]: ph3 = Phono3py(unitcell, supercell_matrix=[3, 3, 2], primitive_matrix='auto', phonon_supercell_matrix=[5, 5, 3])
In [5]: ph3.save("phono3py_disp.yaml")
In [6]: ph3.generate_displacements()Phono3py.generate_fc2_displacements（） 是为 fc2 生成位移数据集的方法。
这通常不需要调用，因为在调用 Phono3py.generate_displacements（） 时会调用它。
如果预期需要对位移模式进行非默认控制，则可以显式调用 Phono3py.generate_fc2_displacements（）。
有关详细信息，请参阅他们的文档字符串。
指定phonon_supercell_matrix时，以下属性可用:1
2
3
4Phono3py.phonon_supercell_matrix
Phono3py.phonon_dataset
Phono3py.phonon_forces
Phono3py.phonon_supercells_with_displacements
·phono3py.load
phono3py.load 的目的是使用从 phono3py-yaml类型文件加载的基本参数创建一个Phono3py类实例，并尝试从当前目录中的 phono3py文件中自动设置力常量和非分析术语校正
在AlN-LDA示例中，单胞结构，超胞矩阵，初级矩阵记录在phono3py_disp_dimfc2.yaml文件中。这很容易借助phono3py.load1
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
50
51
52
53In [1]: import phono3py
In [2]: ph3 = phono3py.load("phono3py_disp_dimfc2.yaml", produce_fc=False)
In [3]: type(ph3)
Out[3]: phono3py.api_phono3py.Phono3py
In [4]: print(ph3.unitcell)
lattice:
- [ 3.110999999491908, 0.000000000000000, 0.000000000000000 ] # a
- [ -1.555499999745954, 2.694205030733368, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 4.978000000000000 ] # c
points:
- symbol: Al # 1
coordinates: [ 0.333333333333333, 0.666666666666667, 0.000948820000000 ]
mass: 26.981539
- symbol: Al # 2
coordinates: [ 0.666666666666667, 0.333333333333333, 0.500948820000000 ]
mass: 26.981539
- symbol: N # 3
coordinates: [ 0.333333333333333, 0.666666666666667, 0.619051180000000 ]
mass: 14.006700
- symbol: N # 4
coordinates: [ 0.666666666666667, 0.333333333333333, 0.119051180000000 ]
mass: 14.006700
In [5]: ph3.supercell_matrix
Out[5]:
array([[3, 0, 0],
[0, 3, 0],
[0, 0, 2]])
In [6]: ph3.nac_params
Out[6]:
{'born': array([[[ 2.51264750e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 2.51264750e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 2.67353000e+00]],
[[ 2.51264750e+00, -2.22044605e-16, 0.00000000e+00],
[ 0.00000000e+00, 2.51264750e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 2.67353000e+00]],
[[-2.51264750e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, -2.51264750e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, -2.67353000e+00]],
[[-2.51264750e+00, 2.22044605e-16, 0.00000000e+00],
[ 0.00000000e+00, -2.51264750e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, -2.67353000e+00]]]),
'factor': 14.39965172592227,
'dielectric': array([[4.435009, 0. , 0. ],
[0. , 4.435009, 0. ],
[0. , 0. , 4.653269]])}
■ Example
@https://github.com/phonopy/phono3py/tree/master/example
AlN-LDA
在这个例子中，我们可以看到较大的fc2超级单元贡献很小，这意味着3x3x2超胞足够好，可以获得良好的声子谱形状。
前提：VASP（LDA+500eV）
超胞的选择：
·fc3 选择3×3×2超胞
·fc2 选择5×5×3超胞
k-mesh:
·原胞：6×6×4
·fc3超胞：2×2×2
·fc2超胞：1×1×2
步骤：
1)计算力
创建完美的具有位移的超胞：phono3py --dim="3 3 2" -c POSCAR-unitcell -d
FORCES_FC3和FORCES_FC2是通过从完美的具有位移的超胞中减去完美超胞的残余力（residual forces）而获得的，即：phono3py --cf3 disp-{00001..01254}/vasprun.xml --cfz disp-00000/vasprun.xml假设 disp3-00000/vasprun.xml 和 disp2-00000/vasprun.xml包含完美超胞的力。在理想情况下，这些力为零，但通常不是。在这里，这被称为"残余力（residual forces）"。
也就是完美的没有位移的超胞原则上是没有力的，但这里包含了，称为残余力（residual forces）,需要减去
目录里的FORCES_FC3被压缩成FORCES_FC3.lzma，解压后可以得到fc3.hdf5 和fc2.hdf5phono3py --sym-fc
2)计算晶格热导
（13x13x9 sampling mesh）
晶格热导可以通过下面来计算：phono3py --mesh="13 13 9" --fc3 --fc2 --br
kappa-m13139.hdf5 被写入为结果。
晶格导热系数计算为k_xx=228.2，k_zz=224.1 W/m-K，300K时
3)非解析校正
使用 —nac 选项，应用非分析项校正，读取 BORN 文件中的 Born 有效电荷和介电常数：phono3py --mesh="13 13 9" --fc3 --fc2 --br --nac
这将300K的热导改变至k_xx=235.7，k_zz=219.1
声子谱的形状对于能量和动量守恒很重要
使用较大的fc2超胞可能会改变声子带结构的形状。
要查看它，请先使用 —dim-fc2 选项重新生成 phono3py_disp.yamlphono3py --dim="3 3 2" --dim-fc2="5 5 3" -c POSCAR-unitcell -d
为了创建力常数，需要FORCES_FC2phono3py --cf2 disp-{00001..00006}/vasprun.xml --cfz disp-00000/vasprun.xml
然后重新创建力常数计算热导1
2phono3py --sym-fc
phono3py --mesh="13 13 9" --fc3 --fc2 --br --nac
Si-cyrstal
关于Si的热导计算
Si-LDA
Si-PBE类似
Si-PBEsol 可以画声子DOS
前提：
超胞是2x2x2的惯用晶胞
VASP：300eV，LDA，不是以Gamma点为中心的2x2x2的k-mesh取样
Si晶体是面心，因此当从惯用晶胞变换到初基晶胞时需要有一个变化矩阵
步骤：
1)产生phono3py_disp.yaml
phono3py -d --dim="2 2 2" -c POSCAR-unitcell --pa="F"
2)创建fc3.hdf5和fc2.hdf5
phono3py --sym-fc
3)利用11x11x11的取样网格来计算晶格热导
phono3py --mesh="11 11 11" --fc3 --fc2 --br
kappa-m111111.hdf5写入结果，计算可以得到：300K的晶格热导112.5 W/m-K
取样是19x19x19，热导是W/m-K
4)累计的热导可以通过phono3py-kaccum脚本求和
phono3py-kaccum --pa="F" -c POSCAR-unitcell kappa-m111111.hdf5 |tee kaccum.dat
为了测试FORCES_FC3的产生：phono3py --cf3 vasprun_xmls/disp-{00001..00111}/vasprun.xml
Si-QE
准备：
QE,超胞是2x2x2 的惯用晶胞
为了获得力：50Ry,2x2x2的k点取样网格
PBE
晶格参数存放在.in文件
Si的晶体是面心结构，因此需要一个变换矩阵ATOMIC_POSITIONS {crystal}一定是要crystal格式的CELL_PARAMETERS {angstrom}一定是要angstrom格式的1
2
3
4CELL_PARAMETERS {angstrom}
3.8564786911 0.0000000000 0.0000000000
-1.9340042258 3.3364748978 0.0000000000
0.0000000000 0.0000000000 17.406650543
步骤：
1)产生phono3py_disp.yaml
phono3py --qe -d --dim="2 2 2" -c Si.in --pa="F"
2)创建fc3.hdf5和fc2.hdf5
phono3py --sym-fc
3)计算热导
11x11x11的网格phono3py --mesh="11 11 11" --fc3 --fc2 --br
kappa-m111111.hdf5写入结果，计算可以得到：300K的晶格热导118.9 W/m-K
19x19x19的网格:300K的晶格热导129.9 W/m-K
4)结果计算
输出结果.out文件都在supercell_out.tar.lzma压缩包里phono3py --cf3 supercell_out/disp-{00001..00111}/Si.out
