あけましておめでとうございます。橋本です。

バンド図は材料物性の動的特性を知る上で重要な情報です。電子はどの面で電気伝導しやすいのか、どの面と接合すれば良いデバイスが作れるか等、特に半導体材料の分野においてはなくてはならないものです。（多分）

最近はARPESで直接バンド図を得る事ができるようになりましたが、（https://arpes.phys.tohoku.ac.jp/contents/study.html ）数値計算による予測は手頃であり、第一原理計算によるシミュレーションはまだまだ現役であるといえます。
さらに、第一原理計算パッケージVASPではグリーン関数を用いて精度を高める計算手法が提供されており、精力的な研究が進められています。

本記事はVASPの出力ファイルである「EIGENVAL」ファイルからバンド図を描くPythonスクリプトを制作しましたので紹介いたします。

#0. 前提
結果の中身について手動で整理する事に興味が無い方は、例えば、pymatgen「https://pymatgen.org/ 」やp4vasp「http://www.p4vasp.at/#/ 」を利用すれば中身をブラックボックスでバンド図を出力する事ができます。

pymatgenについては公式のリファレンスやRKS WEBSITE様「http://ryokbys.web.nitech.ac.jp/pymatgen.html 」の記事が参考になります。
p4vaspはVASP公式サイト「https://www.vasp.at/wiki/index.php/Fcc_Si_bandstructure 」に方法が記載されております。

とは言っても、使い勝手が良いかと言われれば、正直微妙でして。特にバンドが100本を超えるとうまく描けません。本記事では「そんなに難しいフォーマットでは無いので、EIGENVALをpythonで抽出して、手動で描けるようになる」を目標とします。

#1. EIGENVALを得る
バンド図を得るには計算を2回(もしくは3回)走らせる必要があります。

##Step1.1回目の計算
計算に必要なファイルを準備します。

・「POSCAR」：
構造最適化後、又は計算したい構造を作成する等をし、構造が決まったファイルを用意します。

・「KPOINTS」：
計算の精度が収束する十分な数のk点を切ったファイルを用意します。

・「POTCAR」：
好きな擬ポテンシャルファイルを用意します。

・「INCAR」：
計算設定ファイルINCARに用いるそれぞれのタグの意味は「https://qiita.com/youkihashimoto3110/items/de92172e0b5e9f3872d3 」をご参考ください。

1回目の計算は電荷最適化(Electronic Relaxation)計算を新規で回します。構造最適化のパラメータは全てコメントアウトしておきましょう。
必須タグは

ISTART  =   0
PREC    =   Accurate
EDIFF   =   1E-5
LREAL   =   Auto
ALGO    =   VeryFast

です。なお、金属の場合は

ISMEAR  =   1 #または2,-5 
SIGMA   =   0.5

絶縁体・半導体の場合は

ISMEAR  =   0 #または-5
SIGMA   =   0.01

を使用し、outcarのエントロピー項が原子当たり1meV未満になるようにSIGMA値を調整します。
絶縁体・半導体の場合にISMEAR>0を使用しないこと！
また、ISMEAR=-5(テトラヘドロン法)を使用すると精度が向上しますが、大抵の構造の場合適切ではありません。(???)
また、適宜「収束しない場合の有効措置」を使用します。

設定が完了したら1回目の計算を実行します。
##step2. 2回目の計算
計算が終了したら2回目の計算にはいります。2回目の計算では再指定したkメッシュで電荷最適化を継続ジョブで行うため、「INCAR」ファイル、そして「KPOINTS」を書き換えます。

1回目の計算ファイルからの変更点は

ISTART  =   1
ICHARG  =   11
ISYM    =   0
LORBIT  =   11

です。他は同じにして再計算しましょう。
それぞれの意味はテンプレートに記載しています。

「KPOINTS」ファイルは、最悪そのままでも大丈夫ですが、対称性の良いブリュイアンゾーン上でメッシュでを切ったほうが現象が見えやすいと思います。
手動で指定しても良いのですが、ぶっちゃけ2次元だけで勘弁なものを3次元でやる気にはならないので、pymatgenのちからを借ります。
Anacondaを「https://www.anaconda.com/products/individual 」からダウンロード、インストールし、pymatgenをcondaコマンドconda install --channel conda-forge pymatgenでインストールします。
pythonファイルを作成し、以下を実行します。

generate_kpoints.py

from pymatgen.io.vasp.inputs import Kpoints
from pymatgen.core import Structure
from pymatgen.symmetry.bandstructure import HighSymmKpath

struct = Structure.from_file("POSCAR")
kpath = HighSymmKpath(struct)
kpts = Kpoints.automatic_linemode(divisions=40,ibz=kpath)
kpts.write_file("KPOINTS_nsc")

POSCARを読み取り、対称性の良いk点を通るKメッシュ「KPOINT_nsc」が生成されます。
この「KPOINT_nsc」をあたらしい「KPOINTS」を使用します。

準備ができたらジョブを開始しします。

#2. EIGENVALを読み解く
2回目の計算で出力された「EIGENVAL」を任意のテキストエディタで開きます。
EIGENVALの中身は下図ののようになっています。

見えづらくて申し訳ありません。
また、上で示した例は計算でスピンを考慮していません。スピンを考慮した計算の場合、各スピンごとにも出力されます。

#3. pythonで結果を出力する
すなわち、
「固定文字列が6行あり、逆空間の座標3つと？？がある。？？の数値はK点の座標が異なっても変わらない。その下にバンド番号、エネルギー、占有率が次のK座標まで続く。そして次の逆空間の座標3つと？？がある。？？の数値はK点の座標が異なっても変わらない。その下にバンド番号、エネルギー、占有率が...」の構造であるとわかります。また、空改行、空白をスキップする必要があります。

これらを元にプログラミングを行うと以下になります。

band.py

import csv
import matplotlib.pyplot as plt
import numpy as np
import random

eigen_name = "EIGENVAL"           #VASPアウトプットのEIGENVALファイル名
kpoint_name = "KPOINTS"           #VASPのKPOINTSファイル名
save_csv_name = "eigenvalue.csv"  #固有値の保存ファイル
save_band = "band.csv"            #バンド図を保存する

#ファイル名text_nameを指定してリストを返す
def file_to_list(text_name):
    with open(text_name, newline='') as texts:
        read = csv.reader(texts, delimiter=' ', skipinitialspace=True)
        lists = [i for i in read]
    return lists

#保存ファイル名csv_saveとリストlistsを指定してリストを保存
def list_to_csv(csv_name,lists):
    with open(csv_name, mode='w', newline='') as csvname:
        write = csv.writer(csvname)
        write.writerows(lists)

#リストall_listの開始点indexを指定し、グループを切り出し配列ごとに格納
def list_groups(all_list,index):
    groups_all = []
    while True:
        index = index + 1
        temp = []
        count = 0
        temp.append(all_list[index])
        while True:
            index = index + 1
            count = count + 1
            temp.append(all_list[index])
            if len(all_list)==index+1:
                break
            if len(temp[count]) != len(all_list[index+1]):
                groups_all.append(temp)
                index = index + 1 #空白改行を飛ばすため
                break
        if len(all_list)==index+1:
            break
    return groups_all

#kpointを自動でメッシュ切りした時
def k_auto_band(groups_all,k_file):
    spin_up = []
    spin_do = []
    for i in range(len(groups_all)):
        temp_up = []
        temp_do = []
        for j in range(len(groups_all[i])):
            if j == 0: #先頭列にK点を入れる
                k_vec = [float(f) for f in groups_all[i][j][0:3]]
                label = k_file_search(k_file,k_vec)
                for k in groups_all[i][j]:
                    label = label +";"+str(k)
                temp_up.append(label)
                if len(groups_all[i][j])==4:
                    temp_do.append(label)
            else:
                temp_up.append(groups_all[i][j][1])
                if len(groups_all[i][j])==4:
                    temp_do.append(groups_all[i][j][2])
        spin_up.append(temp_up)
        if len(groups_all[i][j])==4:
            spin_do.append(temp_do)
    return spin_up, spin_do

def k_file_search(k_file,lists):
    if k_file[1][0]=="0":
        return ""
    else:
        temp = []
        for i in range(len(k_file)):
            if len(k_file[i]) >= 4:
                for j in range(3):
                    temp.append(float(k_file[i][j]))
                if temp == lists:
                    return k_file[i][4]
                else:
                    return ""


#plotする。Excelはyの系列が256個までしか表示できない。
def graf_show(lists):
    fig = plt.figure()
    list_t = np.array(lists).T
    x = range(len(list_t[0]))
    for i in range(len(lists[0])):
        plt.plot(x,list_t[i])
    fig.savefig("img.png")

#なんか↑のメソッドでも描画できないので200系列をランダムに選ぶ
def random_choice(lists):
    list_t = np.array(lists).T
    random_c = random.sample(range(len(list_t[1:])), 200)
    new_list = []
    new_list.append(list_t[0])
    for i in random_c:
        new_list.append(list_t[i])
    return np.array(new_list).T.tolist()

#メイン
eigen_list = file_to_list(eigen_name)
kpoint_list = file_to_list(kpoint_name)
list_to_csv(save_csv_name,eigen_list)

all_eigen_list = list_groups(eigen_list,6)
spin_up, spin_do = k_auto_band(all_eigen_list,kpoint_list)
list_to_csv(save_band,spin_up)
#graf_show(spin_up)                   #matplotlibでプロット（不具合）
#random_c = random_choice(spin_up)    #200本超えのバンドがある時、バンドからランダムに200本選ぶ
#list_to_csv("radom.csv",random_c)    #↑でランダムに選んだバンド図を保存する.

ひどいプログラムですが、読みにくくは無いと思います。
あとはExcel等でプロットしてあげればバンド図をプロットすることができます。

また、matplotlibモジュールを用いて800バンド程度をプロットしようとしたら処理が止まってしまったので、より軽い方法を模索中です。
Excelでは250系列までしかプロットできないので、200バンドをランダムに抽出してcsv形式で保存するメソッドを用意しました。200バンドを超える系のときrandom_c = random_choice(spin_up)とlist_to_csv("radom.csv",random_c)のコメントアウトを外してください。

上にプログラムを用いた例を示します。