使用GULP(General Utility Lattice Program)计算声子结构和晶格动力学下热导率

GULP(General Utility Lattice Program)基于经典的势能模型，能够以此计算材料的一系列物性，包括晶格优化，原子位置优化，二阶弹性系数，经典势能模型参数拟合，声子计算，缺陷计算等。计算准确程度依赖于势能模型和参数的准确程度，其结果可以作为分子动力学计算的补充。

GULP的官网地址：GULP - Home，其程序安装可参考其manual。

这里主要介绍其利用GULP实现声子计算和晶格动力学下热导率的计算：

GULP计算主要通过输入文本中的keywords和options两大类参数控制，声子计算关键词只需加入phonon即可，热输运的计算有三种形式，一种需要alamode的安装，关键词alamode，alamode可以用于计算原子间高阶力常数，进而计算声子的寿命，并基于运动学声子气的近似得到热导率；一种gulp计算，关键词thermalconductivity，主要用于有缺陷的体系，考虑谐振近似下大量缺陷的局域效应；一种是利用shengBTE，通过gulp计算得到shengBTE所需的CONTROL, FORCE_CONSTANTS_2ND 和FORCE_CONSTANTS_3RD文件：

Calculate the thermal conductivity of a solid at a given temperature. This method uses a gamma point phonon calculation for a supercell to estimate the thermal conductivity using the method of Allen and Feldman, Phys. Rev. B, 48, 12581 (1993). The output includes the mode diffusivities, Di (in cm^2/s), and the overall thermal conductivity (in W/(m.K)). Note that in this approximation the precise value obtained depends on the degree of broadening of the density of states, controlled by the broaden_dos option, the drop tolerance for the Lorentzian broadening function, the temperature used, and finally the size of the supercell (until converged with increasing size). Also note that the non-analytic correction to the phonons at gamma is turned off for a thermal conductivity calculation. In order to correctly allow for the acoustic mode contribution to thermal conductivity then a lower bound frequency should be specified for the Allen-Feldman contribution and an analytic integration used, as described under the omega_af option.
NB: Here the heat capacity for the modes is computed using the full quantised expression from lattice dynamics, rather than approximating the value from equipartition, as used in some other work.
If the thermal conductivity is required more accurately and for ordered solids then GULP can write out the files required for use with the program ShengBTE that implements the Boltzmann
Transport Equations. To obtain the thermal conductivity via this route the steps are:
1) Run a GULP calculation to optimise the structure and compute the gamma point phonons with the option "output shengbte" specified. Note that the third order force constants will be computed analytically if possible, or numerically if they are not available. To force the use of numerical third derivatives the keyword "num3" can be specified. It is important to note that the range of third order force constants output depends on the size of the supercell specified using the "super" option.
2) Modify the file CONTROL to ensure that the correct temperature k point sampling (ngrid) and broadening are specified for your system.
3) Use the files CONTROL, FORCE_CONSTANTS_2ND and FORCE_CONSTANTS_3RD as the inputs for ShengBTE.
4) Assuming the run is successful then the final overall thermal conductivity will be given in the file BTE.kappa_scalar as the last line (converged value from final iteration). Alternatively GULP can use Alamode (if installed) to compute the thermal conductivity and phonon lifetimes. To do this it is important to set the environment variable ALAMODE_DIR to point to the root directory of Alamode.

这里先来介绍gulp基于Alamode的计算：首先需要安装Alamode，这个程序包不依赖平台，将晶胞内原子相互作用的势能泰勒展开，不仅可以基于经典势能模型，也可基于第一性原理计算的结果得到高阶力常数，可以到4阶或5阶。

Alamode官网：ALAMODE — ALAMODE 1.3.0 documentation，其程序独立安装比较费劲，但是可以按照官网推荐的Anaconda中conda结合cmake快速实现安装，相对比较简单。

alamode的主要功能类似于phononpy和phonon3py，但是可以计算4阶，5阶原子间作用的力常数IFCs，以及相应的声子线宽和RTA下的热输运计算，所以程序包含两大模块：alm和anphon，前者用于构建超胞下考虑对称性后有限位移的构型，并利用DFT如QE计算结果给出的力来计算各阶原子间力常数，同phononpy基本功能一致。anphon则用于计算声子色散，态密度和热输运。

在gulp中调用alamode计算经典势能模型下热导率，只需要在gulp的输入文件中：

关键词：optimize(结构优化)  alamode(使用alamode计算热导率)

选项参数中除了加入结构参数，势能模型参数以外，还需要加入：

supercell 5 5  1  #即扩胞，如果上面输入的只是原胞，这个命令可以扩成5*5*1的超胞

alamode用的选项：

ala_disp 0.005 0.01 #即2阶和3阶力常数计算中原子位移量的幅值，注意单位angstrom，与alamode程序默认的Bohr有区别

ala_shrink 16 16 1   #即倒空间k点的密度

temperature 50 50  6  #热导率对应的温度，即起点50K，间隔50K，计算到350K

以上即可，使用gulp运行之前，需要输出alamode安装路径的root目录，如：

export ALAMODE_DIR="alamode安装的根目录"

然后运行gulp，gulp会调用alamode的alm制作超胞的位移pattern，并计算2阶和3阶力常数，并最终用anphon计算RTA近似下上述温度序列下的热导率，计算结果在gulp_rta.results中。

注意：alamode只能实现RTA下热导率的计算，所以意义不是很大，建议直接采用shengBTE或者almaBTE计算收敛后的热导率结果。上述gulp调用alamode的alm，在计算2阶和3阶力常数时，会制作alm的输入文件，内部MODE= fitting，但是alamode在版本1.1.0后将这个模式改为了optimize，所以gulp运行时会出现error：提示不能计算力常数。这时候需要安装alamode的1.0.0版本，即可解决问题。alamode1.0.0版本还没有实现cmake安装，所以安装比较困难，建议只安装alm模块即可，剩余的将高版本alamode的文件拷入即可。

以上由于需要安装旧版本的alamode，因此最好的方法还是自己编程实现接口的文件格式转换，利用gulp作为力计算的引擎，主要流程如下：

1.基于经典势能模型(gulp)建立超胞模型

2.将上述原子结构模型转为qe计算格式文件

3.利用alamode内有处理qe文件的python工具，实现2阶和3阶力常数计算所需的位移pattern和相应的disp***.pw.in系列文件

4.利用python将上述文件格式转为gulp计算格式，gulp计算力需要用到gradient关键词，同时原子坐标后flag需要为1 1 1，计算得到输出文件中力的单位为eV，经过测试需要除相应方向的晶格常数，然后将单位转为Ry/bohr；

5.利用python将gulp计算结果中的力和位移等集中到alamode的DFSET文本文件中，之后就可以用alamode计算声子色散以及RTA下热导率等。

如果需要以上python程序和脚本，可以跟我的qq：1796634232联系索取。