Wannier90とQuantum espressoを使って、タイトバインディングモデルの各ボンド間の重なり積分を見積もることができます。このデータを使って実際にタイトバインディングモデルのバンド図を書いてみます。

Wannier90をつかったタイトバインディングモデル作成は

を参照してください。デモとしてここで紹介されている例でコードを実行します。

バンド計算で確認

出力結果に{prefix}_hr.datというファイルがあり、これに各ボンドでの重なり積分の値が出力されます。

以下ののPythonコードでこれを読み込みバンドを計算することができます。

まずパッケージの読み込み部分です：

calc_bald.py

import numpy as np
import scipy.linalg as LA
import pandas as pd
import time
import matplotlib.pyplot as plt
from matplotlib.path import Path
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 ] )

次に_hr.datファイルの情報をお読み込んでDataframeに直す関数です。

calc_bald.py

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 )

次に入力パラメーターです。

calc_band.py

hr_file = "Bi2Se3_hr.dat"
n_skip = 40
prefix = "Bi2Se3"
output_band = prefix + "_wf_band.csv"

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  ] ])

G = [ 0, 0, 0 ]
Z = [ 0.5, 0.5, 0.5 ]
F = [ 0.5, 0.5, 0 ]
L = [ 0.5, 0, 0 ]
k_symm = np.array([ G,Z,F,G,L])

各パラメーターは以下の通りです。

hr_file：_hrファイルのディレクトリです。
n_skip：Bi2Se3_hr.datの最初の要らない数十行の行数
output_file：バンド図のデータの出力ファイルのディレクトリです。
n_band_id：軌道の総数
fermi：フェルミエネルギー
bravais：基本並進ベクトル(格子長＝1に規格化したもの)、ndarray形で与える。
k_symm：バンドの対称点の列

補足ですが、Bi2Se3_hr.datの最初は以下のようになっていて、重なり積分の値が出力され始めるのは41行目からです。

Bi2Se3_hr.dat

 written on 26May2022 at 17:05:30 
          30
         551
    2    1    1    1    1    1    1    2    4    2    1    2    1    1    2
    1    1    1    2    1    1    1    2    2    1    1    1    1    2    1
    1    1    1    1    2    1    1    1    1    1    1    4    1    1    1
    2    1    2    1    1    1    2    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    2    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    2    1    1    1    1    1    1    2    1    1    1
    2    1    1    2    1    1    1    1    2    1    1    1    1    1    1
    2    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    2    1    1    1
    1    1    2    1    1    1    2    1    2    1    1    1    2    1    1
    1    1    1    2    1    1    1    1    1    1    1    2    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    2
    1    1    1    1    1    1    1    2    1    1    1    1    2    1    1
    2    1    2    2    1    1    1    2    1    1    1    1    1    1    2
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    2    1    1    1    1    1    1    1    2    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    2    1    1    1
    1    1    1    2    1    1    1    2    2    1    2    1    1    2    1
    1    1    1    2    1    1    1    1    1    1    1    2    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    2
    1    1    1    1    1    1    1    2    1    1    1    1    1    2    1
    1    1    2    1    2    1    1    1    2    1    1    1    1    1    2
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    2    1    1    1    1
    1    1    2    1    1    1    1    2    1    1    2    1    1    1    2
    1    1    1    1    1    1    2    1    1    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    2    1    1    1    1    1    1
    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
    1    1    1    1    2    1    1    1    2    1    2    1    1    1    4
    1    1    1    1    1    1    2    1    1    1    1    1    2    1    1
    1    1    2    2    1    1    1    2    1    1    1    2    1    1    2
    1    2    4    2    1    1    1    1    1    1    2
   -5   -1    3    1    1   -0.001467   -0.000000
   -5   -1    3    2    1   -0.000025    0.000043
   -5   -1    3    3    1    0.000995   -0.000257
   -5   -1    3    4    1   -0.000499    0.000038
   -5   -1    3    5    1   -0.001721   -0.000149
   -5   -1    3    6    1    0.000039   -0.000543
   -5   -1    3    7    1    0.003613    0.000000
   -5   -1    3    8    1    0.000000   -0.000000
:
:
:

次にエネルギー準位を計算する波数点の列を用意します。

calc_band.py

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 k_xyz_arr_to_k_crys_arr( k_xyz_arr, bravais ):
    recip = reciprocal( bravais )
    return np.array([ np.dot( LA.inv(recip).T, k ).tolist() for k in k_xyz_arr ])

def add_k_on_way( k_arr0, N ):
    k_list = []
    for i in range( len(k_arr0) - 1 ):
        dk = ( k_arr0[i+1]-k_arr0[i] )/N
        for n in range(N):
            k_list += [ (k_arr0[i] + n*dk).tolist() ]
    k_list += [ k_arr0[-1].tolist() ]
    return np.array(k_list)


def calc_k_trace( k_arr, bravais ):
    recip = reciprocal( bravais )
    k_sum = 0
    k_trace_list = [0]
    for i in range( 1, len(k_arr) ):
        k_sum += LA.norm( np.dot(recip.T,k_arr[i]-k_arr[i-1]), 2 )
        k_trace_list += [k_sum]
    if( len(k_trace_list) > 1 ): k_trace_list = [ K/k_trace_list[-1] for K in k_trace_list ]
    return k_trace_list

重なり積分のDataframeを使ってバンドを計算する関数を作ります。

calc_band.py

def make_hamiltonian_list( k_arr, df_hop, size, fermi ):
    hami_list = [ np.zeros((size,size),dtype=np.complex128) for i in range( len(k_arr) ) ]
    for n_row,sr in df_hop[ abs(df_hop["hop"])>1e-4 ].iterrows():
        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 - fermi*np.identity(size) for h in hami_list ]

def calc_band( k_arr, df_hop, size, fermi, bravais ):
    start = time.time()
    hami_list = make_hamiltonian_list( k_arr, df_hop, size, fermi )
    #band_list = [ LA.eigvalsh( h ).tolist() for h in hami_list ]
    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 ) )
    return pd.DataFrame( band_list, index = calc_k_trace( k_arr, bravais ) )

次のように実行します。

calc_band.py

#hrファイルの読み込み
df_hop = hr_file_to_df_hop_with_band_index( hr_file, n_skip )
# add intermidiate points between symmetric points in the k space
k_arr = add_k_on_way( k_symm, 40 ) #2つ目の引数は対称点の間で補間するk点の個数です。必要に応じて変更してください。今回は40点ずつ補間しました。
#バンド計算
df_band = calc_band( k_arr, df_hop, size, fermi, bravais )
#エネルギーバンドのデータをcsvファイルに書き出し
df_band.to_csv(output_band)

僕の手元のPCでは7並列して12分くらい計算時間がかかりました。このプログラムは(PCのコア数)-1だけ並列処理するように設定してあります。

出力

結果を

で用いたplot_band.pyで第一原理計算の結果と並べて出力してみます。タイトバインディングモデルの結果を出力する処理のコメントアウトを外してください。

(拡大図)

青い点がタイトバインディングモデルから計算したバンドです。完全に一致しているわけではありませんが、概ね合っています。
もしもっと一致させたい場合は、ここから色々とタイトバインディングモデル構築の工夫をしていきます。思いつくものとしては

one-shotではなく、最適化をしてみる。
無視していたBiのs軌道などを入れてみる。
dis_win_min/maxやdis_froz_min/maxを変えてみる。
VESTAなどで原子の結合方向をみて、局所座標を設定してみる。
k点のメッシュをもっと細かくとってみる。

などが挙げられます。
今回はただのデモンストレーションなので、この程度の一致で勘弁してください(汗)。

コメント

タイトバインディングモデルのハミルトニアンを得ることができたので、バンド図でプロットしたエネルギー固有値だけではなく波動関数も取得できます。よって、バンドの混成の具合やスピンの期待値、ベリー位相など様々な量を計算することができるようになります。(ちなみにベリー位相を計算するモードはWannier90に内蔵されています。)