でQuantum Espressoでの電荷密度の計算を紹介しました。ここではそのデータを使って状態密度を各原子軌道の成分に分解したものを計算したいと思います。上記の記事と同様に$\mathrm{Bi_2Se_3}$の例を用います。

状態密度の計算

下準備の計算：

Bi2Se3.nscf.pr.in

&CONTROL
                 calculation = 'nscf' ,
                  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 = 86000,
  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.00000001 ,
		    adaptive_thr = .false. ,
                 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

本計算：

Bi2Se3.pr.in

&PROJWFC
 prefix = 'Bi2Se3'
 outdir = './output'
 filpdos = 'Bi2Se3'
degauss = 0.01
/

filpdosはpdosのデータファイルのprefixです。
degaussは計算結果がスパイキーになるのを防ぐために入れました。ものによっては入れなくても滑らかなものもあります。

これらの計算は以下のように実行します：

../bin/pw.x <Bi2Se3.nscf.pr.in >Bi2Se3.nscf.pr.out
../bin/projwfc.x <Bi2Se3.pr.in >Bi2Se3.pr.out

すると、Bi2Se3.pdos_totとBi2Se3.pdos_atm#...という名前で始まるたくさんのファイルが出力されます。
Bi2Se3.pdos_atm#...はざっくり言うとエネルギー$E$をもつ電子がどこの原子のどの軌道にいるかの割合を表しています。これらの出力結果をまとめてグラフにするコードを作ります。

状態密度の出力

まず、Bi2Se3.pdos_atm#...のファイルをpdosというフォルダに入れておきます。そして次のコードをpdosフォルダの１つ上の階層で実行します。(読み込むファイル名などは適宜変更してください)。

read_pdos.py

import numpy as np
import pandas as pd
import glob
import re

###########################################################################
# move pdos datafile of Quantum Espresso in a folder named "pdos" before! #
###########################################################################

### input ###
prefix = r"Bi2Se3"
Fermi_level = 8.7764

### functions ###
def read_pdos_data(file):
    df = pd.read_csv(file,sep='\s+')
    colname = df.columns.values[1:]
    df = df.dropna(axis=1)
    str2nan = lambda x: np.nan if type(x) == str else x
    df = df.applymap(str2nan).fillna(0)
    #colname[0] = "E(eV)"
    print(file,df.head(),colname,df.columns,"\n")
    df.columns = colname
    return df
    
def exteract_part(name,key):
    return re.search(key, name).group()

def name_to_keyword(filename,dict_key):
    dict1 = {"file":filename}
    for t in dict_key.items():
        flase = exteract_part(filename,t[1])
        dict1[t[0]] = flase
    return dict1

def tag_files(prefix,dict_key):
    filelist = glob.glob(prefix+'**')
    dict_list = []
    for filename in filelist:
        dict_file = name_to_keyword(filename,dict_key)
        dict_list.append(dict_file)
    return dict_list    

def splice_key(dict_f):
    dict_f["orb"] = re.sub(r"wfc#[0-9]+","",dict_f["orb"])
    dict_f["atom"] = re.sub(r"_w","", re.sub(r"atm#[0-9]+","",dict_f["atom"]) ).translate(str.maketrans({'(': '', ')': ''}))
    return dict_f

def get_pdos_df(prefix,EF):
    dict_key = {"orb":r'wfc#\d+(.+)',"atom":"atm#\d+(.*)_w"}
    dict_list = tag_files(prefix,dict_key)
    dict_list = [ splice_key(dict_f) for dict_f in dict_list ]
    dict_df = {}
    for dict_f in dict_list:
        name = dict_f["atom"] + re.sub(r"_j[0-9]+.[0-9]+","",dict_f["orb"])
        print(name)
        #print(dict_f["file"],"\n")
        if( name not in dict_df ): dict_df[name] = read_pdos_data(dict_f["file"])
        else: dict_df[name]["ldos(E)"] += read_pdos_data(dict_f["file"])["ldos(E)"]
    step = 0
    for t in dict_df.items():
        if(step == 0 ):
            df_pdos = pd.DataFrame(t[1]["ldos(E)"])
            df_pdos.columns = [t[0]]
            df_pdos["E"] = t[1]["E(eV)"] - EF
        else: df_pdos[t[0]] = t[1]["ldos(E)"]
        step += 1
    return df_pdos


### execute ###
pdos_file_key = r"pdos/" + prefix + r".pdos_atm#"
df = get_pdos_df(pdos_file_key ,EF = Fermi_level)
df=df.sort_index(axis=1)
df.to_csv("pdos.csv")

すると、各原子軌道の寄与のデータをまとめたものがpdos.csvとして出力されます。
これを次のコードでグラフ描画します。

plot_pdos.py

import pandas as pd
import glob
import re
import os
import matplotlib.pyplot as plt
from matplotlib import colors as clr

df = pd.read_csv("pdos.csv",index_col=0)
df=df.sort_index(axis=1)

plt.rcParams["font.size"] = 13
fig = plt.figure(figsize=(3,5))
ax = fig.add_subplot(111)
for col in df.columns.drop("E"):
    ax = df.plot(x=col,y="E", label=col, ax=ax)
plt.ylim([-2,2])
plt.xlim([0,10])
plt.ylabel("E [eV]")
plt.xlabel("DOS")
plt.axhline(y=0,color= "black",ls="--",lw=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=13)
fig.savefig( "pdos.Bi2Se3.pdf", bbox_inches="tight", pad_inches=0.0, dpi=600)

計算結果がこちらになります。

フェルミ準位の周りはSeのp軌道とBiのp軌道が作るバンドになっていることが分かります。(Biのs軌道も若干混ざってますが、Seのp軌道が支配的なので、まあ、入れなくてもいいでしょう。)

少し物理の話をすると、$\mathrm{Bi_2Se_3}$はラシュバ型スピン軌道相互作用というものが生じる物質として知られています。ラシュバ型スピン軌道相互作用は反転対称性の異なる複数の軌道からなるモデルにおいてスピン軌道相互作用を摂動として取り入れると現れる効果です。この場合のそうした軌道はBiとSeのp軌道の組み合わせが相当します。