Phono3py是一个主要用python写的声子-声子相互作用相关性质的模拟包,可以基于有限位移算法实现三阶力常数和晶格热导率的计算过程,同时输出包括声速,格林奈森常数,声子寿命和累积晶格热导率等参量。
相关介绍和安装请参考往期推荐。
phono3py快速安装教程
QE+phono3py 热导率计算
VASP+phono3py:快速计算晶格热导率
理论到实践:VASP+Phono3py计算Phonon Lifetime
phonopy 和 phono3py 安装教程Phonopy-Spectroscopy计算材料红外和Raman光谱
在phonopy和phono3py的计算过程中,会将过程文件和结果文件写入到hdf5的格式文件中,以此来方便程序读入和写入,phonopy有关的诸多输出量可以通过额外设置获得文本文件,并用于数据处理,phono3py则不能自主输出文本数据。但是在phono3py介绍网址上有如何读取hdf5文件的说明。
http://phonopy.github.io/phono3py/hdf5_howto.html
本文则以此来介绍读取phono3py计算完晶格热导率后的数据读出。
生成有限位移文件后,并通过VASP计算获得力常数fc3.hdf5
和fc2.hdf5后,进行热导率计算,采用11×11×11的网格,
phono3py-load --mesh 11 11 11 --br
有关计算结果会被写入到 kappa-m111111.hdf5中,首先通过pip安装 ipython
pip install ipython
随后打开ipython编辑器
ipython
随后输入(In [1]:是ipython的输入提示,只需要输入冒号后面的部分就行)
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5")
In [3]: list(f)
则输出kappa-m111111.hdf5中包含的数据的title
Out[3]:
['frequency',
'gamma',
'group_velocity',
'gv_by_gv',
'heat_capacity',
'kappa',
'kappa_unit_conversion',
'mesh',
'mode_kappa',
'qpoint',
'temperature',
'weight']
获得热导率张量
In [4]: f['kappa'].shape
Out[4]: (101, 6)
In [5]: f['kappa'][:]
Out[5]:
array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 2.11702476e+05, 2.11702476e+05, 2.11702476e+05,
6.64531043e-13, 6.92618921e-13, -1.34727352e-12],
[ 3.85304024e+04, 3.85304024e+04, 3.85304024e+04,
3.52531412e-13, 3.72706406e-13, -7.07290889e-13],
...,
[ 2.95769356e+01, 2.95769356e+01, 2.95769356e+01,
3.01803322e-16, 3.21661793e-16, -6.05271364e-16],
[ 2.92709650e+01, 2.92709650e+01, 2.92709650e+01,
2.98674274e-16, 3.18330655e-16, -5.98999091e-16],
[ 2.89713297e+01, 2.89713297e+01, 2.89713297e+01,
2.95610215e-16, 3.15068595e-16, -5.92857003e-16]])
获得温度的张量
In [6]: f['temperature'][:]
Out[6]:
array([ 0., 10., 20., 30., 40., 50., 60., 70.,
80., 90., 100., 110., 120., 130., 140., 150.,
160., 170., 180., 190., 200., 210., 220., 230.,
240., 250., 260., 270., 280., 290., 300., 310.,
320., 330., 340., 350., 360., 370., 380., 390.,
400., 410., 420., 430., 440., 450., 460., 470.,
480., 490., 500., 510., 520., 530., 540., 550.,
560., 570., 580., 590., 600., 610., 620., 630.,
640., 650., 660., 670., 680., 690., 700., 710.,
720., 730., 740., 750., 760., 770., 780., 790.,
800., 810., 820., 830., 840., 850., 860., 870.,
880., 890., 900., 910., 920., 930., 940., 950.,
960., 970., 980., 990., 1000.])
获得特定温度热导率和模式热导率(需提前定义)
In [7]: f['kappa'][30]
Out[7]:
array([ 1.09089896e+02, 1.09089896e+02, 1.09089896e+02,
1.12480528e-15, 1.19318349e-15, -2.25126057e-15])
In [8]: f['mode_kappa'][30, :, :, :].sum(axis=0).sum(axis=0) / weight.sum()
Out[8]:
array([ 1.09089896e+02, 1.09089896e+02, 1.09089896e+02,
1.12480528e-15, 1.19318349e-15, -2.25126057e-15])
计算声子寿命
In [9]: g = f['gamma'][30]
In [10]: import numpy as np
In [11]: g = np.where(g > 0, g, -1)
In [12]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
声子群速和热容等数据处理方式同理。
具体修改请查看h5py使用方法和hdf5格式文件读取以及ipython命令使用方法。