3.2 热输运2
三阶力常数矩阵
需要用到thirdorder程序。
使用方法：
- 准备好POSCAR,INCAR,KPOINTS,POTCAR和提交任务脚本（以sub.sh为例）；
其中POSCAR与二阶力常数矩阵求解时的POSCAR-unit一致，KPOINTS也与二阶一致。INCAR文件与二阶力常数求解时的有限位移法一致，其他不变。完成后运行thirdorder实现扩胞。
~/安装目录/thirdorder_vasp.py sow 2 2 2 -12
2 2 2是扩胞倍数，要和二阶力常数矩阵保持一致或者略小。-12指的是考虑多少个近邻原子，需要测试，但是如果计算资源紧缺的时候取个-10就行了，但有风险；正值是指截断半径。之后生成成功以后会出现一堆3RD开头的文件，那些都是POSCAR，一般有几百上千个。
- 参考phonopy计算声子谱的有限位移法，需要把这些POSCAR单独计算。也就是有多少POSCAR，建立多少文件夹，分别进行多少次计算。
find job* -name vasprun.xml|sort -n|/ thirdorder_vasp.py reap a b c -da b c和-d与前面一致。该命令可获得FORCE_CONSTANT_3RD也即是三阶力常数矩阵。
脚本方法
测试-d参数（秦光照老师）
{.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#!/bin/bash
#
# Construct the 3nd files based on the finished calculations
#
# Written BY QIN, GuangZhao <qin.phys@gmail.com>
# 2016-11-22
DIR="/vol-th/home/xhaiyan/WangN/Au6S2/third241-14"
if [ ! -d "$DIR" ];then
echo Please specify the source directory.
exit
fi
for NAME in 3RD.POSCAR.*
do
NUM="${NAME#3RD.POSCAR.}" # the number of the 3nd POSCAR file
CODE=$(head -1 $NAME) # the code of the POSCAR file
# find the source POSCAR file in the source directory
INFO=$(grep $CODE $DIR/3RD.POSCAR.* | head -1)
[ x"$INFO" == x ] && echo "NOT FOUND for job-$NUM" && continue
SOURCE_NUM="${INFO%:*}" # the number of the source POSCAR file
SOURCE_NUM="${SOURCE_NUM##*.}"
echo "$DIR/job-$SOURCE_NUM/ --> job-$NUM/"
# diff 3RD.POSCAR.$NUM $DIR/3RD.POSCAR.$SOURCE_NUM
rsync -a $DIR/job-$SOURCE_NUM/vasprun.xml ./job-$NUM/
done针对大量POSCAR批量建立文件夹
{.line-numbers} 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15#module load spglib/master
#thirdorder_vasp.py sow 3 3 1 -10
#生成2*2*2的超胞 0.6表示截断半径为6埃。
#sow a b c m abc表示生成的超胞大小，m为正数的时候表示截断半径（单位：nm）（例如0.6，表示截断半径为6.0埃）；为负数时表示第几近邻截断（此时应取负整数，例如-3表示第三近邻）
for i in 3RD.POSCAR.*
do
s=$(echo $i| cut -d "." -f 3)
d=job-$s
mkdir $d
cp $i $d/POSCAR
cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d
pushd $d
popd
done批量提交任务
{.line-numbers} 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17#!/bin/sh
root_path='pwd'
#在每个job文件夹里进行静态计算。
for i in {001..1050} #按顺序从job-001到1050做静态vasp计算,根据生成的文件夹设定
do
cd job-$i
echo "job-$i", time
sbatch sub.sh
wait
cd ..
done
#上述静态计算完成后提取三阶力常数文件
#cd ${root_path}
#find job* -name vasprun.xml| sort -n| python thirdorder_vasp.py reap 2 2 1 -12 #最后四个数字跟你扩胞时一致
脚本赋予权限：chmod +x .sh;运行脚本：./sh
ShengBTE获得热输运性质
运行ShengBTE计算热输运性质
- 获得二阶力常数计算：
计算流程与声子谱计算流程相同，用DFPT和有限位移法都可以。将得到的FORCE_CONSTANTS重命名为FORCE_CONSTANTS_2ND。 - 获得三阶力常数计算：见上，输出为FORCE_CONSTANTS_3RD文件。
- 准备CONTROL文件：格式如下:
{.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#CONTROL for ShengBTE
&allocations
nelements=2 #元素种类
natoms=6 #原子个数
ngrid(:)=20 20 5 #计算网格密度
&end
&crystal
lfactor=0.1 #软件长度单位换算
lattvec(:,1)=3.1613337552503711 0.0000000000000025 0.0000000000000000 #原胞晶胞参数
lattvec(:,2)=-1.5806668776251651 2.7377953418477841 0.0000000000000000
lattvec(:,3)=-0.0000000000000000 0.0000000000000000 12.3318234981499426
elements= "Mo" "S" #原子顺序
types= 1 1 2 2 2 2 #原子个数，2代表第二种元素，4个2代表有4个第二类原子
positions(:,1)= 0.3333333449999998 0.6666666890000030 0.2500000000000000 #原子坐标，按顺序
positions(:,2)=0.6666666240000012 0.3333333229999980 0.7500000189999980
positions(:,3)= 0.3333333449999998 0.6666666890000030 0.6230719397502427
positions(:,4)=0.6666666240000012 0.3333333229999980 0.3769280602499874
positions(:,5)=0.6666666240000012 0.3333333229999980 0.1230720067501656
positions(:,6)=0.3333333449999998 0.6666666890000030 0.8769280222497613
scell(:)=4 4 1 #扩胞倍数，和前面一致
&end
¶meters
T=300 #温度
&end
&flags
nonanalytic=.TRUE.
nanowires=.FALSE.
convergence=.TRUE. #迭代法，默认最大1000步
&end
详细参数可参考手册。
- 提交任务
- 获得不同温度的数据，只需要依次修改CONTROL文件的T tag。或者通过脚本：
1 | for i in {350..1350..100} # 对i赋值，区间为300到900，每搁10取一个数，即（300 310 320 330 ....890 900），可以根据自身情况自行设置 |
数据处理
- 计算完成后，先检查收敛情况，找到T300K文件夹下的BTE.kappa,第一列代表迭代步数，一般未达到默认1000步可认为收敛正常。
- 计算完成后在当前目录下会生成T300K文件夹和一些输出文件。文件处理参见天玑算b站视频：https://www.bilibili.com/video/BV1Yh411n7tx/?spm_id_from=333.788&vd_source=17084cc867ff3bb86398ff1a79b443f2
其他参考：https://www.bilibili.com/read/cv22090733/
注意事项
计算得到的热导率和文献报道有较大偏差，可能偏大或者偏小，可能的原因有：
1.采用的2阶力常数FORCE_CONSTANT中超胞内原子排序不是按照ucatomsc_zsc_y*sc_x排列；力常数文件与control文件中原胞内原子次序是否一致；
2.shengBTE的control文件中晶格参数单位默认nm，有转换系数的参数；
3.使用QE利用phonopy计算得到的FORCE_CONSTANT文件中，单位需要转换，不能直接利用shengBTE计算，vasp利用phonopy则单位正确；
3阶力常数的计算量较大，但精度要求很高，如0.0003Rh/bohr^2量级的力的误差，也会导致最终热导率的数值出现大的偏离；如使用cp2k计算力常数，通常使用QE等平面波软件，对于超胞体系，随着扩胞，在较大的截断半径下，所需计算量非常大，使用CP2K中的QS模块也能进行第一性原理计算，具有相对较高的计算效率，但是CP2K不同于平面波的QE等，计算中仅有gamma点计算，需要采用实空间较大超胞，如果使用和QE等相同的超胞，则会导致原子受力精度不高，一般在~0.0003Rh/bohr^2量级的误差，也会对某些体系，特别是二维体系的热导率有很大的影响。在使用cp2k进行3阶力常数大量计算之前，先进行超胞收敛测试，可以在相同的位移构型下，不断扩胞，判断力计算的变化情况和误差。用QE等计算也建议先进行k点等计算参数的收敛判断，否则会浪费大量时间，热导率计算结果却有很大偏离。
使用经验势能模型，在通过lammps进行热输运计算之前，可以先利用分子静力学软件GULP结合shengBTE模块，预判势能应用于热输运计算的准确性。但经验势能用于热输运精度和前景一般，机器势能表现更好。
原文链接：https://blog.csdn.net/odin_linux/article/details/130084321
静态计算可以各个job文件夹同时分别计算，不用按顺序来。等所有job文件夹计算完后，自行运行“find job…” 提取三阶力常数的命令。严格来说，胞的大小和三阶力常数计算的截断半径都应该进行测试，测试以所计算得到的热导率的收敛进行判断。计算资源有限无法测试的话，尽可能取大一点。Born有效电荷和介电常数的计算应该不是必要的，如果你要考虑在Gamma点附近LO-TO劈裂的话(主要针对离子晶体)，则需要进行计算。可以加上，感觉要好点。
二阶力常数和三阶力常数所用的超胞大小可以不统一，我看文献中有这种情形，不过一般作者采用的二阶力常数的超胞都是大于等于三阶力常数超胞的大小。
原文链接：https://zhuanlan.zhihu.com/p/147810223
四声子过程
除了ShengBTE的CONTROL文件中的常规输入，FourPhonon还需要四阶力常数和一些新的标签在CONTROL文件中。请检查ShengBTE网站以获取其他标签的定义。
并行环境
1.1和1.0版本是为MPI并行性编写的。从1.2版本开始，支持四声子散射的迭代求解器，我们已经迁移到OpenMP，以处理该迭代求解器所需的大内存。未来我们可能会支持MPI+OpenMP混合并行性，以允许更多灵活性。确保在编译时添加-qopenmp，或者：
export FFLAGS=-qopenmp -traceback -debug -O2 -static_intel
四阶IFCs文件：FORCE_CONSTANTS_4TH
该文件包含四阶原子间力常数矩阵，并使用稀疏描述来节省空间。要构建四阶力常数，可以参考Fourthorder的python脚本。该力常数的格式是对三阶力常数的直接扩展，包含nb个这样的矩阵块：
- 一行空白
- 一个从1开始的顺序索引（从1到nb）
- 一行表示第二个晶胞的笛卡尔坐标（单位：Å）
- 一行表示第三个晶胞的笛卡尔坐标（单位：Å）
- 一行表示第四个晶胞的笛卡尔坐标（单位：Å）
- 一行表示四个相关原子的从1开始的索引，每个原子索引从1到natoms
- 81行的力常数矩阵（单位：eV Å⁴）。开头的索引用于标记笛卡尔坐标轴。
该文件的一个示例：
1 |
|
CONTROL 文件设置
此 CONTROL 文件包含所有用户指定的设置和参数，包括晶体结构信息、展宽因子、q 网格、温度、函数等。
一般用法
要启用 FourPhonon 功能，您需要添加一个新的 &flags 名称列表：
1 | four_phonon （逻辑，默认为 .false.）：计算四声子相空间和四声子散射率。 |
标志的实际用法
以下是这些标志的一些实用组合：
1 | onlyharmonic=.true. 和 four_phonon=.true.：仅计算四声子相空间。 |
注意： CONTROL 文件中的所有其他参数，例如温度或 q 网格，适用于三声子和四声子过程。FourPhonon 包不支持纳米线功能。
加速 RTA 解的采样方法
为了从散射过程样本中估计散射率，请引用 [Z. Guo 等人，npj Comput. Mater. 10, 31 (2024).] 并使用相关标签。
参数设置
num_sample_process_3ph 和 num_sample_process_3ph_phase_space
num_sample_process_3ph 和 num_sample_process_3ph_phase_space （整数，默认值为 -1）：表示从每个模式中提取的样本数，以用于估计三声子相空间和散射率。 -1 表示不进行取样，执行严格计算。注意，当 convergence=.true.，即使用迭代方案计算三声子散射时，num_sample_process_3ph 必须为 -1，因为取样方法基于弛豫时间近似。
num_sample_process_4ph 和 num_sample_process_4ph_phase_space
num_sample_process_4ph 和 num_sample_process_4ph_phase_space （整数，默认值为 -1）：表示从每个模式中提取的样本数，以用于估计四声子相空间和散射率。 -1 表示不进行取样，执行严格计算。注意，当 four_phonon=.false. 时，num_sample_process_4ph 和 num_sample_process_4ph_phase_space 必须均为 -1。
示例
计算三声子散射率（精确计算，迭代），并在 RTA 层面采用 100,000 个样本过程计算四声子散射：
1 | ¶meters |
TDEP IFCs 的输入
对于 TDEP 格式的谐波力常数，请使用以下描述的标签/脚本：
tdep=.true. 启用 TDEP 格式谐波 IFCs 的处理，该输入文件名应始终为 infile.forceconstant。如果系统是极性材料，除了常规参数（如 Born 有效电荷和介电常数）外，我们还需要用户手动输入耦合参数 (\lambda_{Ewald})，该参数在 TDEP IFCs 的末尾近似写出，如下所示：
1 | 0.675614929199219 # Coupling parameter in Ewald summation |
注意：输入的谐波 IFCs 是未带极性修正的裸 IFCs，以便 FourPhonon 能够处理并在内部添加修正。该模块在很大程度上依赖 TDEP 源代码，我们无法保证成功。建议用户检查处理后的频率和三声子散射率。
对于非极性材料，我们还开发了一种脚本，将 TDEP 格式转换为 Phonopy 格式，名为 tdep2phonopy.f90。这样您无需激活 tdep 标志。编译 Fortran 后，运行可执行文件（以下是示例）：
1 | ./tdep2phonopy outfile.forceconstant 10 10 1 2 |
这将 TDEP IFCs （outfile.forceconstant）转换为 Phonopy 格式，构建一个来自双原子单位胞的 (10 \times 10 \times 1) 超胞。请检查转换后的声子色散。我们建议使用足够大的超胞来存储所有的第二阶 IFC 中的成对相互作用。
示例
计算三声子散射率（精确计算，迭代），并在 RTA 层面计算四声子散射。使用的 TDEP 格式谐波 IFCs 并应用非解析修正，Ewald 参数为 (\lambda =0.675614929199219)。同样需要 Born 有效电荷（此处未列出），与 ShengBTE 对极性系统的处理相似。
1 | ¶meters |
TDEP 格式的 3阶/4阶 IFCs
对于 TDEP 格式的 3阶/4阶 IFCs，我们开发了一个独立的脚本，将其转换为 ShengBTE 格式，名为 tdep2ShengBTE.f90。我们需要在同一目录下有 infile.ucposcar。编译 Fortran 后，运行可执行文件：
示例
1 | ./tdep2ShengBTE 3 2 outfile.forceconstant_thirdorder |
这将 TDEP 的 3 阶 IFCs（outfile.forceconstant_thirdorder）转换为 ShengBTE 格式。
Output files
Besides the routine output files from previous ShengBTE program, FourPhonon generates these output files:
- BTE.Numprocess_4ph: number of allowed four-phonon scattering processes, for each irreducible q point and phonon band
- BTE.P4: phase space available for four-phonon processes, for each irreducible q point and phonon band
- BTE.P4_total: total volume in phase space available for four-phonon processes
- BTE.P4_plusplus*, BTE.P4_plusminus*, BTE.P4_minusminus*: similar to BTE.P4 but only includes contributions from ++/+-/- - processes
*Note for four-phonon scatterings, there are three different channels: recombination (++), redistribution (+-) and splitting (- -) processes.
Under temperature-dependent directories: (the unit of scattering rates is ps−1 if not specified)
- BTE.WP4: weighted phase space available for four-phonon processes ([rad/ps]$^{-5}$, 2nd column) vs angular frequency (rad/ps, 1st column) for each irreducible q point and phonon band
- BTE.WP4_plusplus*, BTE.WP4_plusminus*, BTE.WP4_minusminus*: similar to BTE.WP4 but only includes contributions from ++/+-/- - processes
- BTE.w_3ph: three-phonon scattering rates for each irreducible q point and phonon band, this file replaces the original output BTE.w_anharmonic. Absorption and emission processes are written out into BTE.w_3ph_plus and BTE.w_3ph_minus
- BTE.w_4ph: four-phonon scattering rates for each irreducible q point and phonon band. Similarly, we provide the contributions from different channels: BTE.w_4ph_plusplus, BTE.w_4ph_plusminus and BTE.w_4ph_minusminus
- BTE.w_4ph_normal: four-phonon scattering rates of normal processes, for each irreducible q point and phonon band. This file has 5 columns and each one represents: angular frequency in rad/ps, recombination channel, redistribution channel, splitting channel, overall scattering rates from normal processes
- BTE.w_4ph_Umklapp: four-phonon scattering rates of Umklapp processes, for each irreducible q point and phonon band. The format is the same as BTE.w_4ph_normal
- BTE.kappa*: thermal conductivity results are written out as usual