使用alamode计算非谐热导

赵文博Zhao

已于 2023-09-05 21:25:17 修改

阅读量1k

点赞数 3

文章标签： python

于 2023-02-14 16:04:24 首次发布

本文链接：https://blog.csdn.net/lielie12138/article/details/129028637

版权

0.软件安装

a.按照官网说的配置conda环境

% conda create --name alamode -c conda-forge python=3
% conda activate alamode
% conda install -c conda-forge compilers openmpi boost eigen cmake spglib fftw scipy numpy h5py ipython
% pip list #检查安装情况
#下载软件
% git clone https://github.com/ttadano/alamode.git
% cd alamode
% git checkout develop
#安装
% mkdir _build
# cmake是一个自动生成Makefile的软件。
% cmake -S . -B _build #-S CMakeLists.txt目录，-B 构建目录
% cmake --build "./"

计算流程

用alm获得位移模式
计算位移过的配置的原子力
通过拟合估计力常数
计算声子色散和声子DOS
估计用于热导的非谐IFCs
计算热导率
分析结果

本文依照此教程进行，各种输入文件在 example/Si 目录

https://alamode.readthedocs.io/en/latest/tutorial.html

1.用alm获得位移模式

alm软件拟合IFC（interatomic force constants）通过最小二乘拟合 **（by least square fitting）**输入文件si_alm.in。

&general
  PREFIX = si222
  MODE = suggest
  NAT = 64; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1  # 1: harmonic, 2: cubic, ..
/

&cell
  20.406 # factor in Bohr unit
  1.0 0.0 0.0 # a1
  0.0 1.0 0.0 # a2
  0.0 0.0 1.0 # a3
/

&cutoff 
  Si-Si None
/

&position
  1 0.0000000000000000 0.0000000000000000 0.0000000000000000   
  1 0.0000000000000000 0.0000000000000000 0.5000000000000000
  1 0.0000000000000000 0.2500000000000000 0.2500000000000000
  1 0.0000000000000000 0.2500000000000000 0.7500000000000000
  1 0.0000000000000000 0.5000000000000000 0.0000000000000000
  1 0.0000000000000000 0.5000000000000000 0.5000000000000000
  1 0.0000000000000000 0.7500000000000000 0.2500000000000000
  ...
  后面的没写哈，例子文件里有完整的坐标信息。一共是64个，这个要和NAT参数保持一致，不然会报错。

执行

alm si_alm.in > si_alm.log1

生成文件si222.pattern_HARMONIC这里是si的建议位移模式。

Basis : C
    1:    1
      1              1              0              0

在模式文件中，建议的位移模式在笛卡尔坐标中定义。如您在文件中看到的那样，单质SI的谐波IFC只有一个位移模式。

2.计算位移过的配置的原子力

接下来，按照SI22222.Pattern_Harmonic提供的位移模式，计算所有位移过的配置的原子力。

为此，您首先需要确定位移的大小ΔU，这应该很小，以使非谐的贡献可以忽略不计。在大多数情况下，ΔU〜0.01Å是一个合理的选择。
然后，准备为每种配置运行外部DFT代码所需的输入文件。由于此过程有些复杂，因此我们为VASP，Quantum-Espresso（QE）和XTAPP提供了一个Python脚本。在tools/ directory中使用脚本sumplace.py，您可以按以下方式生成必要的输入文件：
（这里只给了QE的，实际也可以用VASP等其他软件）

python displace.py --QE=si222.pw.in --mag=0.01 -pf si222.pattern_HARMONIC

注意！要自己准备si222.pw.in这个文件！他是原始的QE输入文件！这个脚本是根据这个原始的输入文件改的。！
K_POINTS automatic
4 4 4 1 1 1
111是必须的！！！不然算不出来位移后的结果！！
结果会生成disp1.pw.in文件，这个就是位移后的输入文件。有多少位移模式会生成多少这个文件。
disp1.pw.in和si222.pw.in文件目前看只会移动原子位置，不会影响别的。

对所有位移后配置执行pw.x

#!/bin/bash

# Assume we have 20 displaced configurations for QE [disp01.pw.in,..., disp20.pw.in].

for ((i=1;i<=20;i++))
do
    num=`echo $i | awk '{printf("%02d",$1)}'`
    mkdir ${num}
    cd ${num}
    cp ../disp${num}.pw.in .
    pw.x < disp${num}.pw.in > disp${num}.pw.out
    cd ../
done

下一步是通过python脚本提取物（也在tools/目录中）收集位移数据和强制数据。该脚本可以从多个输出文件中提取原子位移，原子力和总能量：如下：
(这一步只需要.out文件生成的其他文件都不需要)
(这一步运行如果报错大概是因为numpy的版本太高了，降到1.21就行了)

python extract.py --QE=si222.pw.in *.pw.out > DFSET_harmonic

所有输出文件会被整理到DFSET_harmonic中。数据格式为DFSET。关于这个格式的详细内容可以参考

https://alamode.readthedocs.io/en/latest/almdir/inputalm.html#label-format-dfset

在这里插入图片描述

3.通过拟合估计力常数

编辑文件si_alm.in以执行最小二乘拟合。
更改为MODE = optimize

&general
  PREFIX = si222
  MODE = optimize   # <-- here
  NAT = 64; NKD = 1
  KD = Si
/

增加一个部分

&optimize
  DFSET = DFSET_harmonic
/

再次执行alm

alm si_alm.in > si_alm.log2

在这个过程中alm读取DFSET_harmonic文件，生成si222.fcs，si_alm.log2，si222.xml三个文件。他们的内容分别是：

- si222.fcs

包含所有IFCs在Rydberg原子单位。在这个文件的第一部分我们可以找到对称不可约IFCs集合（symmetrically irreducible sets of IFCs ）

*********************** Force Constants (FCs) ***********************
 *        Force constants are printed in Rydberg atomic units.       *
 *        FC2: Ry/a0^2     FC3: Ry/a0^3     FC4: Ry/a0^4   etc.      *
 *        FC?: Ry/a0^?     a0 = Bohr radius                          *
 *                                                                   *
 *        The value shown in the last column is the distance         *
 *        between the most distant atomic pairs.                     *
 *********************************************************************

 ----------------------------------------------------------------------
      Index              FCs         P        Pairs     Distance [Bohr]
 (Global, Local)              (Multiplicity)                           
 ----------------------------------------------------------------------

  *FC2
       1       1     2.7613553e-01   1     1x     1x       0.000
       2       2    -9.5252308e-04   2     1x     2x      10.203
       3       3    -3.7574247e-04   2     1x    33x      10.203
       4       4     8.8561868e-03   1     1x     3x       7.215
       5       5     1.9633030e-03   1     1x     3y       7.215
       6       6    -3.9798276e-03   1     1x    17x       7.215
       7       7    -3.9080855e-03   1     1x    17z       7.215
       8       8     8.2188178e-03   4     1x     6x      14.429
       9       9     4.3878651e-05   4     1x    34x      14.429
      10      10    -6.7358757e-02   1     1x     9x       4.418
      11      11    -4.7360832e-02   1     1x     9y       4.418
      12      12     1.0670028e-03   1     1x    10x       8.460
      13      13    -8.2098635e-04   1     1x    10y       8.460
      14      14    -6.9742830e-04   1     1x    10z       8.460
      15      15     2.1135669e-04   1     1x    26x       8.460
      16      16    -4.3031422e-03   1     1x    11x      11.118
      17      17     5.0759340e-04   1     1x    11y      11.118
      18      18    -4.7178438e-04   1     1x    25x      11.118
      19      19    -4.7518256e-04   1     1x    25z      11.118
      20      20    -7.4138963e-04   2     1x    20x      12.496
      21      21     1.2900836e-03   2     1x    20y      12.496
      22      22    -1.7263732e-04   2     1x    20z      12.496
      23      23     6.4515756e-05   2     1x    35x      12.496
      24      24     3.5369980e-04   1     1x    28x      13.254
      25      25    -5.1486930e-04   1     1x    28y      13.254

第三列：超胞的简谐 IFCs Φμν(i,j)

第四列：the multiplicity P，他是每个相互作用（I，J）发生在给定的截止半径内的次数。例如对于第二行pairs(1x,2x)，P=2，这是因为r1,2=r1,2′ ，其中2`是周期性边界条件下原子2的相邻图像（完全没懂）

如果比较IFC的大小，则第三列中的值应该除以P

- si_alm.log2

包含fitting error

https://alamode.readthedocs.io/en/latest/almdir/formalism_alm.html#fitting-formalism

在这里插入图片描述

-si222.xml

包含晶体结构，对称性，IFC和随后的声子计算所需的所有其他信息。

4.计算声子色散和声子DOS

-计算声子

si_phband.in

&general
  PREFIX = si222
  MODE   = phonons
  FCSXML = si222.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

&cell
  10.203
  0.0 0.5 0.5
  0.5 0.0 0.5
  0.5 0.5 0.0
/

&kpoint
  1  # KPMODE = 1: line mode
  G 0.0 0.0 0.0 X 0.5 0.5 0.0 51
  X 0.5 0.5 1.0 G 0.0 0.0 0.0 51
  G 0.0 0.0 0.0 L 0.5 0.5 0.5 51
/

参数FCSXML需要指定上一步计算的.xml文件。
参数&cell中不再是超胞的晶格参数，而是primitive cell的晶格参数。
在声子分散计算中，＆kpoint 中的第一个条目应为1（kpmode = 1）。

执行

anphon si_phband.in > si_phband.log

生成文件si222.bands，这个文件里是延路径的声子频率。单位:cm-1.
你可以用gnuplot或任何其他绘图软件绘制声子色散关系。为了可视化声子色散关系，我们在tools/目录中提供了Python脚本plotband.py（需要Matplotlib）

python plotband.py si222.bands

您可以在此窗口中将图形保存为png、eps或其他格式。还可以通过–unit选项将声子频率的能量单位从cm-1更改为THz或meV。有关plotband.py用法的更多详细信息，请键入

python plotband.py -h

不过，这个脚本算的很慢，直接把si222.bands文件拿出来用origin画会快很多
在这里插入图片描述

在这里插入图片描述
两个画的一模一样

-计算DOS

复制si_phband.in为si_phdos.in，更改&kpoint部分

&kpoint
  2  # KPMODE = 2: uniform mesh mode
  20 20 20
/

执行：

anphon si_phdos.in > si_phdos.log

生成文件：

si222.dos：phonon DOS文件
si222.thermo：thermodynamic functions热力学函数

为了可视化声子DOS和投影的DOS，我们在tools/目录中提供了Python脚本plotdos.py（需要Matplotlib）

python plotdos.py --emax 550 --nokey si222.dos

要提高DOS的分辨率，请使用更密集的k网格和更小的DELTA_E值重试。
在这里插入图片描述
这里，如果是没有图形输出功能的机器可以去plotdos.py里面把plt.show()这行改成plt.savefig(‘dos.png’)。这样就不会输出图形，而是保存成dos.png文件。

5. 估计用于热导的非谐IFCs

将文件si_alm.in复制到si_alm2.in。按如下方式编辑si_alm2.in的&general、&interaction和&cutoff字段：

&general
  PREFIX = si222_cubic
  MODE = suggest
  NAT = 64; NKD = 1
  KD = Si
/

&interaction
  NORDER = 2
/

将NORDER标记从NORDER=1更改为NORDER=2，以计算立方IFC（这里立方指三次方，指考虑三阶力常数）。在这里，我们通过将截止半径指定为7.3来考虑到次近邻的立方相互作用对

&cutoff
  Si-Si None 7.3
/

这个数字是怎么确定的呢？
7.3bohr 略大于第二近邻的距离（7.21461bohr）。可以根据自己的情况适当更改截止值。（原子距离可以在文件sialm.log中找到。）（为啥要大于第二临近距离？）
在这里插入图片描述
执行

alm si_alm2.in > si_alm2.log

生成文件
si222_cubic.pattern_HARMONIC：还是简谐的可能位移与第一步获得的si222.pattern_HARMONIC文件一样。

si222_cubic.pattern_ANHARM3：考虑三阶非谐的可能位移模式

然后，计算文件si222_cubic.pattern_ANHARM3中给出的位移配置的原子力，并将位移和力数据集收集到文件DFSET_cubic中，就像在步骤3和4中对谐波IFC所做的那样。（这一步省略了，就是把这几种模式都拿去自洽，最后会生成一个DFSET_cubic文件。）

在ai_alm2.in中
更改MODE = suggest 为 MODE = optimize，并增加下面部分。

&optimize
  DFSET = DFSET_cubic
  FC2XML = si222.xml # Fix harmonic IFCs
/

FC2XML标记简谐IFC所在文件-si222.xml

执行

alm si_alm2.in > si_alm2.log2

生成文件
si222_cubic.fcs
si222_cubic.xml.
内容和之前算的简谐的是一样的，参见目录3.1和3.3，不过这里是算了三阶的。

注意：
前面我们是固定的简谐力常数的情况下计算了包含三阶力常数的力常数（约等于在二阶的基础上加三阶）。我们也可以不固定简谐力常数。直接生成三阶力常数。

把二阶和三阶的dft的位移结果合并

cat DFSET_harmonic DFSET_cubic > DFSET_merged

将读取的文件更改为合并的文件

&optimize
  DFSET = DFSET_merged
/

6.计算热导率

在 si_phdos.in的基础上改输入文件 si_RTA.in

&general
  PREFIX = si222
  MODE = RTA
  FCSXML = si222_cubic.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

换更粗一些的k网格

&kpoint
  2
  10 10 10
/

执行

anphon si_RTA.in > si_RTA.log

请耐心等待。这可能需要一段时间。作业完成后，您可以找到一个文件si222.kl，其中保存了晶格热导率。您可以使用gnuplot（或任何其他绘图软件）绘制此文件，如下所示：

$ gnuplot
gnuplot> set logscale xy
gnuplot> set xlabel "Temperature (K)"
gnuplot> set ylabel "Lattice thermal conductivity (W/mK)"
gnuplot> set term pngcairo  #保存为PNG格式
gnuplot> set output "a.png" #保存到文件
gnuplot> plot "si222.kl" usi 1:2 w lp

也可以直接用origin画图，画前两排，x和y都取log。
在这里插入图片描述

最后画完就是紫色那一条。可以看出紫色和绿色在T趋近于0的地方有分歧。（绿色是准确的，考虑了声子边界散射的图）
这是因为我们在目前的计算中只考虑了本征声子-声子散射，而忽略了在低温范围内占主导地位的声子边界散射。可以使用工具目录中的python脚本analyze_phonons.py来包括边界散射的效果：

cp ~/soft/alamode/tools/analyze_phonons.py ~/soft/alamode/_build/tools
#第一次执行记得把脚本复制到安装目录
python ~/soft/alamode/_build/tools/analyze_phonons.py --calc kappa_boundary --size 1.0e+6 si222.result > si222_boundary_1mm.kl

脚本的运行原理：
在这个脚本中，声子寿命是使用Matthiessen法则改变的**（the Matthiessen’s rule）**
在这里插入图片描述

这里，方程右侧的第一项是声子-声子散射引起的散射速率，第二项是尺寸为L的晶界引起的散射率。必须使用–size选项以nm为单位指定尺寸L。结果也显示在上图中，边界效应解决了T趋近于0时热导率发散的问题。
注意！当使用 smearing method（ISMEAR=0或1）而不是四面体法（ the tetrahedron method ）（ISMEAR=-1）进行计算时，即使没有边界效应，热导率也可能在非常低温的区域具有峰值结构。由于 smearing method中使用的有限拖尾宽度（ the finite smearing width ϵ），因此出现了该峰值。当我们降低ε值时，κ的峰值应消失。此外，非常密集的q网格对于描述低温区域中的声子激发和热输运是必要的（无论ISMEAR值如何）。

7. 分析结果

有许多方法可以分析结果，以更好地理解纳米尺度的热输运。下面介绍了一些特殊的函数：

a.声子寿命

文件si222.result包含不可约k点处的声子线宽。您可以从该文件中提取声子寿命，如下所示：

$ analyze_phonons.py --calc tau --temp 300 si222.result > tau300K_10.dat
$ gnuplot
gnuplot> set xrange [1:]
gnuplot> set logscale y
gnuplot> set xlabel "Phonon frequency (cm^{-1})"
gnuplot> set ylabel "Phonon lifetime (ps)"
gnuplot> set term pngcairo
gnuplot> set output "tau.png"
gnuplot> plot "tau300K_10.dat" using 3:4 w p

在这里插入图片描述

在上图中，用20×20×20q点计算的声子寿命也画出来了。

b.Cumulative thermal conductivity(累计热导率)

定义在这里

https://alamode.readthedocs.io/en/latest/anphondir/formalism_anphon.html#cumulative-kappa

在这里插入图片描述

$ analyze_phonons.py --calc cumulative --temp 300 --length 10000:5 si222.result > cumulative_300K_10.dat
$ gnuplot
gnuplot> set logscale x
gnuplot> set xlabel "L (nm)"
gnuplot> set ylabel "Cumulative kappa (W/mK)"
gnuplot> set term pngcairo
gnuplot> set output "Cum.png"
gnuplot> plot "cumulative_300K_10.dat" using 1:2 w lp

在这里插入图片描述
图中也画出了202020的结果。

c.Thermal conductivity spectrum(热导率光谱)

更改 si_RTA.in文件

&general
  PREFIX = si222
  MODE = RTA
  FCSXML = si222_cubic.xml

  NKD = 1; KD = Si
  MASS = 28.0855

  EMIN = 0; EMAX = 550; DELTA_E = 1.0 # <-- frequency range
/

&cell
10.203
 0.0 0.5 0.5
 0.5 0.0 0.5
 0.5 0.5 0.0
/

&kpoint
 2
 10 10 10
/

&analysis
 KAPPA_SPEC = 1 # compute spectrum of kappa
/

使用EMIN、EMAX和DELTA_E标记指定频率范围，并在&analysis字段中设置KAPPA_SPEC=1。然后，再次执行anphon

anphon si_RTA.in > si_RTA2.log

计算完成后，您可以找到文件si222.kl_spec，其中包含各种温度下的热导率谱。您可以按以下方式绘制室温下的数据：

$ awk '{if ($1 == 300.0) print $0}' si222.kl_spec > si222_300K_10.kl_spec
$ gnuplot
gnuplot> set xlabel "Frequency (cm^{-1})"
gnuplot> set ylabel "Spectrum of kappa (W/mK/cm^{-1})"
gnuplot> plot "si222_300K_10.kl_spec" using 2:3 w l lt 2 lw 2