Quantum Espressoでのバンド計算

https://qiita.com/hiro949/items/2ea86ef94df2f5d7ab41

でQuantum Espressoで電荷密度分布を求めたので、それを使ってバンドを計算したいと思います。同様に$Bi_2Se_3$の例を用いて計算します。

# バンド計算

対称線上のバンド構造を計算します。対称線の取り方は次の論文が参考になると思います：

https://www.sciencedirect.com/science/article/pii/S0927025610002697?via%3Dihub

$\mathrm{Bi_2Se_3}$の対称性は$R\bar{3}m$で、この論文のA.11.のパターンの結晶になっています（$\alpha=60^\circ$）。次のページを参照してください。

http://pd.chem.ucl.ac.uk/pdnn/symm1/rhomb.htm

この対称点は論文中のTable14で与えられます。これを結ぶ線(論文のFig.14)のうち$\Gamma-Z-F-\Gamma-L$を考えます。（これは対称線ではありませんが、この物質は$k_x-k_y$平面に射影して考えることが多いのでこの線分が定番になっています。）

まずは下準備の計算を行います。次のファイルを入力とします。

```Bi2Se3.nscf.band.in
 &CONTROL
                 calculation = 'bands' ,
                  wf_collect = .true. ,
                      outdir = './output' ,
                  pseudo_dir = '~/q-e-qe-6.5/pseudo/' ,
                       prefix = 'Bi2Se2' ,
		      disk_io = 'low',
                   verbosity = 'high' ,
                     tstress = .true. ,
                     tprnfor = .true. ,
  max_seconds = 170000,
  restart_mode = from_scratch,
 /
 &SYSTEM
  ibrav = 0,
  celldm(1) = 7.81213,
  nat = 5,
  ntyp = 2,
                       ecutwfc = 150 ,
                     ecutrho = 1500 ,
                 occupations = 'smearing' ,
                     degauss = 0.00001 ,
		     noncolin = .true. ,
		     lspinorb = .true. , ,
		        nbnd = 240
 
 /
 &ELECTRONS
            electron_maxstep = 10000000,
                    conv_thr = 0.0000000001 ,
                 mixing_beta = 0.5 ,
/
CELL_PARAMETERS {alat}
  0.577350269189626   0.000000000000000   2.306079664570230 
 -0.288675134594813   0.500000000000000   2.306079664570230 
 -0.288675134594813  -0.500000000000000   2.306079664570230 
ATOMIC_SPECIES
  Bi  208.98000  Bi.rel-pbe-dn-kjpaw_psl.1.0.0.UPF
  Se   78.96000  Se.pbe-dn-kjpaw_psl.1.0.0.UPF
  
ATOMIC_POSITIONS {crystal}
Bi   0.400460000000000   0.400460000000000   0.400460000000000 
Bi   0.599540000000000   0.599540000000000   0.599540000000000 
Se   0.209700000000000   0.209700000000000   0.209700000000000 
Se   0.790300000000000   0.790300000000000   0.790300000000000 
Se   0.000000000000000   0.000000000000000   0.000000000000000 


K_POINTS crystal_b
5
0 0 0 20
0.5 0.5 0.5 20
0.5 0.5 0 20
0 0 0 20
0.5 0 0 20
```

Bi2Se3.scf.inとの違いとしては以下の通りです。

- &control-calculationがbandsになっている。
- &systemにnbndが追加されている。これは各波数点で低い順に何個のエネルギー順位を計算するかのパラメーターである。
- K_POINTが計算する波数点の線分の対称点の個数と位置になっている。(逆格子ベクトルの係数と線分間で計算する波数点の個数を入力する。今回は20点ずつサンプリングした。)

その後にバンド計算を行います。その時の入力ファイルは次の通りです。

```Bi2Se3.band.in
&bands
 outdir = './output',
 prefix = 'Bi2Se3',
 filband = 'Bi2Se3.band',
 lsym=.false.,
/
```
filbandは出力ファイルの名前(拡張子抜き)です。
lsymは各対称点で小群の規約表現の指標を求めるモードです。今回は単にバンド図を描きたいだけなのでオフにします。

これらの計算は以下のコマンドで実行します。(もちろん並列化した方がよいです。)

```
../bin/pw.x <Bi2Se3.nscf.band.in >Bi2Se3.nscf.band.out
../bin/bands.x <Bi2Se3.band.in >Bi2Se3.band.out

```

すると、Bi2Se3.band.gnuというgnuplotの画像データが出力されます。gnuplotで出力してもいいのですが、次のコードを使ってpythonで出力します。

### バンド図の出力

```plot_band.py
import pandas as pd

def band_gnu_to_df(gnu_file):
    df = pd.read_csv(gnu_file, delim_whitespace=True, names=('k','E'))
    k_arr = df["k"].unique()
    g = df.groupby("k")
    n_row = len(g.get_group(k_arr[0])["E"])
    df_out =  pd.DataFrame(index=[],columns=list(range(n_row)))
    for k in k_arr:
        df_out.loc[k] = g.get_group(k)["E"].values
    return df_out

#２回の計算に分けた場合にそれぞれの出力されたgnuファイルの結果をつなげる
def merge_band_gnu(gnu_file1,gnu_file2):
    df_band1 = band_gnu_to_df(gnu_file1)
    df_band2 = band_gnu_to_df(gnu_file2)
    df_band2 = df_band2.drop(df_band2.index[0],axis=0)
    df_band2.index = df_band2.index.values + df_band1.index.values[-1]
    df_band = df_band1.append(df_band2)
    return df_band

# 入力パラメーター
gnu_file = "Bi2Se3.band.gnu"
Fermi_level = 8.7764

df_band = band_gnu_to_df(gnu_file) - Fermi_level

# 波数点を規格化
df_band.index = df_band.index.values/df_band.index.values[-1]

# グラフ描画のためのパラメーター
yrange = [-2,2]
symm_point_names = ["Γ","Z","F","Γ","L"]
symm_point_arr = df_band.index.values[::20]
image_file = "band.pdf"

#比較のためのタイトバインディング模型のバンド(後述)
#df_wann = pd.read_csv("Bi2Se3_wf_band.csv",index_col=0)

#グラフ描画
import matplotlib.pyplot as plt
import matplotlib.lines as mlines

fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)

ax = df_band.plot(legend=False,ylim=yrange,color="red",lw=2.,ax=ax)
#ax = df_wann.plot(legend=False,ylim=yrange,color="blue",lw=0,marker="x",ms=2.,ax=ax)

red_line = mlines.Line2D([], [], color='red',label='SCF')
plt.legend(handles=[red_line],fontsize=15, bbox_to_anchor=(1.01, 1), loc='upper left')
#blue_x = mlines.Line2D([], [], color='blue',linestyle=":", label='MLWF')
#plt.legend(handles=[red_line,blue_x],fontsize=15, bbox_to_anchor=(1.01, 1), loc='upper left')

# print symmetric point as x-axis lavels
ax.set_xlabel("")
ax.set_ylabel("Energy (eV)",fontsize=18)
ax.set_xticks(symm_point_arr)
ax.set_xticklabels(symm_point_names)
ax.tick_params(labelsize=18)
ax.axhline(y=0,color='black',linestyle=':')
for k in symm_point_arr:
    ax.axvline(x=k,color='black',linestyle=':')
ax.set_xlim([symm_point_arr[0],symm_point_arr[-1]])

ax.tick_params(direction='in', length=6, width=2,pad=10)
plt.tight_layout()
fig.show()
#fig.savefig( image_file,bbox_inches="tight", pad_inches=0.0, dpi=300)
```


計算結果はこんな感じです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/b3d0d11c-af52-8dcd-470a-de4c61618355.png)

実験でも$0.3$ eVくらいのバンドギャップが空いていると分かっていて、それとも大体合っています。