FermiSurferでフェルミ面を描画しよう

物性物理学ではしばしは材料のに応答特性をイメージするためにフェルミ面の図が必要とされることがあります。
そこでFermisurferというソフトを使ってこのようなプロットを行う方法をご紹介します。

https://fermisurfer.osdn.jp/index.html.ja

今回は第一原理計算の結果から得たFermi面をプロットしたいと思います。

例として$Bi_2Se_3$という物質の計算結果を用いたいと思います。この材料の第一原理計算はこちらの記事をご覧ください。

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

これにはDOS計算(成分分解しないもの)の下処理の計算を行う必要があるので、まずはQuantum espressoでのDOS計算の仕方をお伝えしたいと思います。

# Quantum Espressoでの計算

## DOSの計算の入力

次の入力ファイルを用意します。

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

でのバンド計算の下準備の計算とほとんど同じです。
k点メッシュを$8\times8\times8$にしましたが、フェルミ面を書くにはかなり制度が必要なので、もっと細かく取る必要があるかもしれません。

```Bi2Se3.nscf.dos.in
  &CONTROL
                 calculation = 'bands' ,
                  wf_collect = .true. ,
                      outdir = './output' ,
                  pseudo_dir = '~/q-e-qe-6.5/pseudo/' ,
                       prefix = 'Bi2Se3' ,
		      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 automatic
 8 8 8 0 0 0
```

次にDOS計算の入力ファイルです。

```Bi2Se3.dos.in
&DOS
 outdir = './output',
 prefix = 'Bi2Se3',
 fildos = 'Bi2Se3.dos',
/
```

## フェルミ面計算のための入力

フェルミ面計算の入力は次のようになります。

```Bi2Se3.fs.in
&fermi
 outdir = './output',
 prefix = 'Bi2Se3',
/
```

## 実行

これらを次のコマンドで実行します。

```
../bin/pw.x <Bi2Se3.nscf.dos.in >Bi2Se3.nscf.dos.out
../bin/dos.x <Bi2Se3.dos.in >Bi2Se3.dos.out
../bin/fs.x <Bi2Se3.fs.in >Bi2Se3.fs.out
```

２つ目のコマンドが終了したらBi2Se3.dosというファイルが生成されます。ここにDOSのデータが入っています。

```Bi2Se3.dos
#  E (eV)   dos(E)     Int dos(E) EFermi =    4.367 eV
 -38.638  0.6845E-01  0.6845E-03
 -38.628  0.1377E-80  0.6845E-03
 -38.618  0.1377E-80  0.6845E-03
 -38.608  0.1377E-80  0.6845E-03
 -38.598  0.1377E-80  0.6845E-03
 -38.588  0.1377E-80  0.6845E-03
 -38.578  0.1377E-80  0.6845E-03
 -38.568  0.1377E-80  0.6845E-03
 -38.558  0.1377E-80  0.6845E-03
 -38.548  0.1377E-80  0.6845E-03
 -38.538  0.1377E-80  0.6845E-03
 -38.528  0.1377E-80  0.6845E-03
 -38.518  0.1377E-80  0.6845E-03
 -38.508  0.1377E-80  0.6845E-03
 -38.498  0.1377E-80  0.6845E-03
 -38.488  0.7043E-66  0.6845E-03
 -38.478  0.1377E-80  0.6845E-03
 -38.468  0.1377E-80  0.6845E-03
:
:
:
```

また、3つ目のコマンドが終了すると、Bi2Se3_fs.bxsfというファイルが出力されます。これがフェルミ面のデータになっています。このデータをFermiSurferで読み込んでFermi面を描画したいと思います。


# FermiSurferによるFermi面の描画

次にFermiSurferを使って得られたFermi面のデータを実際に出力してみます。

## インストール

公式ページは

http://fermisurfer.osdn.jp/index.html.ja

です。また以下のページからダウンロードします。

https://ja.osdn.net/projects/fermisurfer/releases/

## 描画

以下、Windows環境で考えます。FermiSurfer.exeをダウンロードします。ただ、これ単独では実行することはできません。Bi2Se3_fs.bxsfをプログラムから開くなどからFermiSurfer.exeで開きます。

すると、次のような図を得ることができます。図はマウスのドラッグやRotateテキストボックスで回転させることができます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/071b551d-78fe-fc41-ec79-9b4d95fe0c29.png)

## 断面

任意の方向で断面を描くこともできます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/02273cc8-e2d7-22e0-2101-e3410819f12e.png)

まず、section-v(黄色の網掛け)で断面の法線の係数(基底は逆格子ベクトル)を入力します。
次にsectionボックス(青色の網掛け)にチェックを入れます。そしてupdate(赤矢印)をクリックすると、上図のように断面を出力できます。この図は$(1,1,1)$方向の図です。
また、sectionボックスの下のSection fileというボタンをクリックすると、断面のデータ(fermi_line.dat：フェルミ面の断面のデータとbz_line.dat：ブリリュアンゾーンの断面の頂点のデータ)が出力されます。

# タイトバインディング模型から求める場合

次のコードでfrmsfファイルを作成してFermiSurferで出力しても、フェルミ面を描画することができます。これは

https://qiita.com/hiro949/items/b9fd1d345b463e298c99

でWannier90を使って求めたタイトバインディング模型で第１ブリリュアンゾーンのエネルギーバンドを求めてそのデータを使っています。ただ、計算にはかなり時間がかかるので気を付けてください。

```calc_frmsf.py
import numpy as np
import scipy.linalg as LA
import pandas as pd
import time

# system size
Lx=8
Ly=8
Lz=8

# band index number
nbnd = 94

# input parameters
output_file = "BiSe3.frmsf"
hr_file = "Bi2Se3_hr.dat"
n_skip = 40
n_band_id = 30
fermi = 8.7764
bravais = np.array([ [   0.577350269189626,   0.000000000000000,   2.306079664570230], 
 [-0.288675134594813,   0.500000000000000,   2.306079664570230,], 
 [-0.288675134594813,  -0.500000000000000,   2.306079664570230  ] ])

# %%
from joblib import Parallel, delayed

def parallel( iteration, Njob, Nver, func  ):
    return Parallel(n_jobs=Njob,verbose=Nver)( [delayed(func)(i) for i in iteration ] )

# %%
def reciprocal( bravais_T ):
    reciprocal = np.zeros((3,3))
    deno = LA.det( bravais_T.T )
    for i in [0,1,2]:
        reciprocal[i] = np.cross( bravais_T[(i+1)%3], bravais_T[(i+2)%3] )/deno
    return reciprocal

# %%
def hr_file_to_df_hop_with_band_index( hr_file, n_skip ):
    col_name_list = ["a1","a2","a3","n_band1","n_band2","hop_real","hop_imag"]
    col_dict = {"a1":int,"a2":int,"a3":int,"n_band1":int,"n_band2":int,"hop_real":float,"hop_imag":float}
    df_hop = pd.read_table(hr_file, skiprows=[ i for i in range(0,n_skip)], header=None,
                       sep="\s+", names= col_name_list, dtype=col_dict)
    df_hop["hop"] = df_hop["hop_real"] + 1j*df_hop["hop_imag"]
    df_hop["n_band1"] -= 1; df_hop["n_band2"] -= 1
    return df_hop.drop( ["hop_real","hop_imag"], axis=1 )

def make_hamiltonian_list( k_arr, df_hop, size, fermi, threshold=1e-4 ):
    hami_list = [ np.zeros((size,size),dtype=np.complex128) for i in range( len(k_arr) ) ]
    E0 = fermi*np.identity(size)
    for sr in df_hop[ abs(df_hop["hop"])>threshold ].itertuples():
        an = np.array([ sr.a1, sr.a2, sr.a3 ])
        t = sr.hop
        row = sr.n_band1
        col = sr.n_band2
        for i, k in enumerate( k_arr ):
            phase = 2*np.pi*np.dot(an,k)
            hami_list[i][row,col] += t*( np.cos(phase) + 1j*np.sin(phase) )        
    return [ h - E0 for h in hami_list ]

def calc_band( Lx, Ly, Lz, df_hop, size, fermi ):
    start = time.time()
    n_arr = np.array([ [i1,i2,i3] for i1 in range(Lx) for i2 in range(Ly) for i3 in range(Lz)])
    k_arr = np.array([ [n[0]/Lx,n[1]/Ly,n[2]/Lz] for n in n_arr ])
    hami_list = make_hamiltonian_list( k_arr, df_hop, size, fermi )
    band_list = Parallel(n_jobs=-2,verbose=4)( [delayed(LA.eigvalsh)(h) for h in hami_list ] )
    band_list = [ Ek.tolist() for Ek in band_list ]
    end = time.time()
    print( "{0} sec.".format( end - start ) )
    
    df = pd.DataFrame(columns=['x','y','z','n','Energy'])
    df_band = pd.DataFrame(n_arr,columns=['x','y','z']).join(pd.DataFrame(band_list))
    for i in range(size):
        df0 = df_band.loc[:,['x','y','z',i]].rename(columns={i:'Energy'})
        df0['n'] = [ i for j in range(Lx*Ly*Lz) ]
        df = df.append(df0)
    return df

# %%
df_hop = hr_file_to_df_hop_with_band_index( hr_file, n_skip )

# %%
df_band = calc_band(Lx, Ly, Lz, df_hop, n_band_id, fermi)

# %%
df_band = df_band.sort_values(['n','x','y','z'])
recip = reciprocal( bravais ).round(6)

# %%
with open(output_file, 'w') as f_handle:
    f_handle.write(str(Lx)+" "+str(Ly)+" "+str(Lz)+"\n")
    f_handle.write("1\n")
    f_handle.write(str(nbnd)+"\n")
    np.savetxt(f_handle, recip, delimiter=' ', fmt ='%.6f')
    np.savetxt(f_handle, df_band["Energy"].values, delimiter=' ', fmt ='%.6f')
```

# 参考
- http://fermisurfer.osdn.jp/index.html.ja
- https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.105.076802
- https://iopscience.iop.org/article/10.1088/1367-2630/12/6/065013