import numpy as np import os import re import sys # sysモジュールを追加 def parse_outcar_fermi_energy(outcar_filepath="OUTCAR"): """ OUTCARファイルからフェルミ準位 (E_fermi) を取得する。 Args: outcar_filepath (str): OUTCARファイルのパス。デフォルトはカレントディレクトリの"OUTCAR"。 Returns: float or None: 取得したフェルミ準位 (eV)。見つからない場合はNone。 """ if not os.path.exists(outcar_filepath): print(f"エラー: '{outcar_filepath}' が見つかりません。") return None e_fermi = None try: with open(outcar_filepath, 'r') as f: # ファイルを逆順に読み込み、最後に出現するE-fermiを探すのがより正確 # ただし、ファイル全体を読み込む必要があるため、大きなファイルでは注意 # ここでは、行を順に読み込み、最後に見つかったものを採用 for line in f: if "E-fermi" in line: match = re.search(r"E-fermi\s*:\s*(-?\d+\.\d+)", line) if match: e_fermi = float(match.group(1)) if e_fermi is None: print(f"警告: '{outcar_filepath}' からフェルミ準位が見つかりませんでした。") return e_fermi except Exception as e: print(f"エラー: '{outcar_filepath}' の解析中に問題が発生しました: {e}") return None def analyze_vasp_bands_from_eigenval(eigenval_filepath="EIGENVAL", outcar_filepath="OUTCAR", custom_fermi_energy=None): """ VASPのEIGENVALファイルとOUTCARファイルからCBM, VBM, およびバンドギャップを取得する。 占有率情報に基づいてバンド端を特定する。金属の場合にはバンドギャップを0とする。 Args: eigenval_filepath (str): EIGENVALファイルのパス。デフォルトはカレントディレクトリの"EIGENVAL"。 outcar_filepath (str): OUTCARファイルのパス。デフォルトはカレントディレクトリの"OUTCAR"。 custom_fermi_energy (float, optional): 手動で指定するフェルミ準位。指定しない場合はOUTCARから取得。 Returns: dict: CBM, VBM, バンドギャップの値、および関連情報を含む辞書。 取得できなかった場合はNoneを返す。 """ if not os.path.exists(eigenval_filepath): print(f"エラー: '{eigenval_filepath}' が見つかりません。") return None # フェルミ準位の取得 e_fermi = None if custom_fermi_energy is not None: e_fermi = custom_fermi_energy print(f"カスタムフェルミ準位を使用: {e_fermi:.3f} eV") else: e_fermi = parse_outcar_fermi_energy(outcar_filepath) if e_fermi is None: print("エラー: フェルミ準位が取得できませんでした。バンドギャップの計算をスキップします。") return None print(f"OUTCARから取得したフェルミ準位 (E_fermi): {e_fermi:.3f} eV") try: with open(eigenval_filepath, 'r') as f: lines = f.readlines() # EIGENVALヘッダーの解析 # 5行目からNIONS, NKPTS, NBANDSを読み取る # EIGENVALの形式: # 0: コメント # 1: 原子数 # 2: 空行 # 3: 空行 # 4: 空行 # 5: NIONS NKPTS NBANDS (VASP 5.x/6.x) # 6: 空行 # 以降、k点ごとにデータブロック # k点座標と重み # NBANDS行: バンド番号 エネルギー 占有率 (NSPIN=1) # バンド番号 エネルギー_up 占有率_up エネルギー_down 占有率_down (NSPIN=2) # 空行 (k点間の区切り) if len(lines) < 6: print("エラー: EIGENVALファイルの行数が少なすぎます。ヘッダーを読み取れません。") return None # NIONS, NKPTS, NBANDSは通常6行目 (インデックス5) にある header_line_5 = lines[5].strip().split() if len(header_line_5) < 3: print("エラー: EIGENVALのヘッダー形式が予期と異なります (5行目)。NIONS, NKPTS, NBANDSが見つかりません。") return None # NIONSは不要だが、形式合わせのため読み込む # nions = int(header_line_5[0]) nkpts = int(header_line_5[1]) nbands = int(header_line_5[2]) # 全てのバンドデータを格納するリスト # 各要素は辞書: {"k_idx", "band_idx", "spin_idx", "energy", "occupancy"} band_data_all = [] nspin = 0 # スピンチャンネル数を初期化、データから決定する is_metallic = False # 金属であるかどうかのフラグ # データ開始行 (0-indexed)。ヘッダーの次の空行 (インデックス6) から開始 current_line_idx = 6 for ik in range(nkpts): # k点間の空行をスキップ (もしあれば) # current_line_idx がファイルの終端を超えないようにチェック if current_line_idx < len(lines) and not lines[current_line_idx].strip(): current_line_idx += 1 # k点座標と重みの行をスキップ # current_line_idx がファイルの終端を超えないようにチェック if current_line_idx >= len(lines): print(f"エラー: EIGENVALファイルの終端に達しました (k点 {ik+1} のデータ不足)。k点座標行が見つかりません。") break current_line_idx += 1 # バンドデータ開始行のチェック if current_line_idx >= len(lines): print(f"エラー: EIGENVALファイルの終端に達しました (k点 {ik+1} のバンドデータ不足)。") break for iband in range(nbands): # 各バンド行の前にファイルの終端チェック if current_line_idx >= len(lines): print(f"エラー: EIGENVALファイルの終端に達しました (k点 {ik+1}, バンド {iband+1} のデータ不足)。") break band_line = lines[current_line_idx].strip().split() current_line_idx += 1 if not band_line: # 空行の場合、スキップまたはエラー print(f"警告: EIGENVALのバンドデータ行が空です (k点 {ik+1}, バンド {iband+1})。スキップします。") continue # 最初のk点、最初のバンドで列数をチェックし、NSPINを決定 if nspin == 0: # nspinがまだ決定されていない場合のみ実行 # band_line[0] はバンド番号 # 残りの列の数でスピンチャンネルを判断 num_data_cols = len(band_line) - 1 if num_data_cols == 2: # エネルギー 占有率 (NSPIN=1) nspin = 1 elif num_data_cols == 4: # エネルギー1 占有率1 エネルギー2 占有率2 (NSPIN=2) nspin = 2 else: print(f"エラー: EIGENVALのバンドデータ行の列数が予期と異なります: {len(band_line)} (バンド番号を含む)。3列または5列を期待しました。") return None print(f"検出されたスピンチャンネル数 (NSPIN): {nspin}") # エネルギーと占有率の抽出 # band_line[0] はバンド番号 energy_up = float(band_line[1]) occupancy_up = float(band_line[2]) band_data_all.append({ "k_idx": ik, "band_idx": iband, "spin_idx": 0, # スピンアップ (0-indexed) "energy": energy_up, "occupancy": occupancy_up }) # 部分占有のチェック (占有率が厳密に0でも1でもない場合、金属と判断) # 微小な数値誤差を考慮して閾値を設定 if abs(occupancy_up - 0.0) > 1e-4 and abs(occupancy_up - 1.0) > 1e-4: is_metallic = True # print(f"デバッグ: 部分占有バンド発見 (k={ik+1}, band={iband+1}, spin=1, occ={occupancy_up:.6f})") # デバッグ用 # 金属と判断されたら、これ以上占有率をチェックする必要はない # ただし、全てのバンドデータを読み込む必要があるので、ループは続ける if nspin == 2: energy_down = float(band_line[3]) occupancy_down = float(band_line[4]) band_data_all.append({ "k_idx": ik, "band_idx": iband, "spin_idx": 1, # スピンダウン (0-indexed) "energy": energy_down, "occupancy": occupancy_down }) # 部分占有のチェック (スピンダウン) if abs(occupancy_down - 0.0) > 1e-4 and abs(occupancy_down - 1.0) > 1e-4: is_metallic = True # print(f"デバッグ: 部分占有バンド発見 (k={ik+1}, band={iband+1}, spin=2, occ={occupancy_down:.6f})") # デバッグ用 # 各k点での全バンドデータを読み込んだ後、current_line_idxは次のk点ブロックの開始、 # またはファイル終端を指している。次のループで空行スキップが処理される。 if not band_data_all: print("エラー: バンドエネルギーのデータが読み込めませんでした。EIGENVALファイルが空であるか、形式が不正です。") return None # バンドギャップの計算 vbm_info = {"energy": -np.inf, "k_idx": -1, "band_idx": -1, "spin_idx": -1} cbm_info = {"energy": np.inf, "k_idx": -1, "band_idx": -1, "spin_idx": -1} # 占有率の閾値 (VASPは通常0か1だが、数値誤差を考慮) # 占有されているとみなす閾値 (例: 0.5より大きい) OCCUPIED_THRESHOLD = 0.5 # 占有されていないとみなす閾値 (例: 0.5より小さい) UNOCCUPIED_THRESHOLD = 0.5 found_occupied_band = False found_unoccupied_band = False for data in band_data_all: energy = data["energy"] occupancy = data["occupancy"] k_idx = data["k_idx"] band_idx = data["band_idx"] spin_idx = data["spin_idx"] # VBMの候補: 占有されているバンド (occupancy > OCCUPIED_THRESHOLD) の中で最大のエネルギー if occupancy > OCCUPIED_THRESHOLD: found_occupied_band = True if energy > vbm_info["energy"]: vbm_info["energy"] = energy vbm_info["k_idx"] = k_idx + 1 # 1-indexed vbm_info["band_idx"] = band_idx + 1 # 1-indexed vbm_info["spin_idx"] = spin_idx + 1 # 1-indexed # CBMの候補: 占有されていないバンド (occupancy < UNOCCUPIED_THRESHOLD) の中で最小のエネルギー if occupancy < UNOCCUPIED_THRESHOLD: found_unoccupied_band = True if energy < cbm_info["energy"]: cbm_info["energy"] = energy cbm_info["k_idx"] = k_idx + 1 # 1-indexed cbm_info["band_idx"] = band_idx + 1 # 1-indexed cbm_info["spin_idx"] = spin_idx + 1 # 1-indexed # VBM/CBMが特定できなかった場合の最終チェック if not found_occupied_band: print("エラー: 占有されたバンドが見つかりませんでした。VBMを特定できません。") return None if not found_unoccupied_band: print("エラー: 占有されていないバンドが見つかりませんでした。CBMを特定できません。") return None vbm_energy = vbm_info["energy"] cbm_energy = cbm_info["energy"] band_gap = 0.0 # 金属判定ロジック if is_metallic: # 部分占有のバンドが見つかった場合 band_gap = 0.0 print("情報: 部分占有のバンドが見つかりました。この系は金属的であると判断し、バンドギャップは0とします。") elif cbm_energy > vbm_energy: # 部分占有がなく、CBM > VBMの場合（絶縁体/半導体） band_gap = cbm_energy - vbm_energy else: # 部分占有がなく、CBM <= VBMの場合（バンドが重なっている金属/半金属） band_gap = 0.0 print("情報: CBMがVBM以下です。この系は金属的または半金属的であると判断し、バンドギャップは0とします。") results = { "NSPIN": nspin, "NKPTS": nkpts, "NBANDS": nbands, "E_fermi": e_fermi, # フェルミ準位は参考情報として含める "VBM_Energy": vbm_energy, "VBM_K_Index": vbm_info["k_idx"], "VBM_Band_Index": vbm_info["band_idx"], "VBM_Spin_Index": vbm_info["spin_idx"], "CBM_Energy": cbm_energy, "CBM_K_Index": cbm_info["k_idx"], "CBM_Band_Index": cbm_info["band_idx"], "CBM_Spin_Index": cbm_info["spin_idx"], "Band_Gap": band_gap } return results except Exception as e: print(f"エラー: '{eigenval_filepath}' の解析中に問題が発生しました: {e}") # 詳細なエラー情報も表示 import traceback traceback.print_exc() return None if __name__ == "__main__": # コマンドライン引数の解析 eigenval_path = "EIGENVAL" outcar_path = "OUTCAR" custom_fermi = None # 引数の順序は EIGENVAL_path, OUTCAR_path, custom_fermi_energy の順を想定 if len(sys.argv) > 1: eigenval_path = sys.argv[1] if len(sys.argv) > 2: outcar_path = sys.argv[2] if len(sys.argv) > 3: # ここを修正しました try: custom_fermi = float(sys.argv[3]) except ValueError: print(f"警告: 無効なフェルミ準位の引数 '{sys.argv[3]}'。OUTCARから取得します。") print(f"'{eigenval_path}' と '{outcar_path}' からバンド情報を解析しています...") band_info = analyze_vasp_bands_from_eigenval(eigenval_path, outcar_path, custom_fermi) if band_info: print("\n--- バンド計算結果 ---") print(f"スピンチャンネル数 (NSPIN): {band_info['NSPIN']}") print(f"k点数 (NKPTS): {band_info['NKPTS']}") print(f"バンド数 (NBANDS): {band_info['NBANDS']}") print(f"フェルミ準位 (E_fermi): {band_info['E_fermi']:.3f} eV") print("-" * 30) print(f"VBM (Valence Band Maximum): {band_info['VBM_Energy']:.3f} eV") print(f" K点番号: {band_info['VBM_K_Index']}") print(f" バンド番号: {band_info['VBM_Band_Index']}") print(f" スピンチャンネル: {band_info['VBM_Spin_Index']}") print("-" * 30) print(f"CBM (Conduction Band Minimum): {band_info['CBM_Energy']:.3f} eV") print(f" K点番号: {band_info['CBM_K_Index']}") print(f" バンド番号: {band_info['CBM_Band_Index']}") print(f" スピンチャンネル: {band_info['CBM_Spin_Index']}") print("-" * 30) print(f"バンドギャップ (Band Gap): {band_info['Band_Gap']:.3f} eV") else: print("\nバンド情報の取得に失敗しました。") print("EIGENVALまたはOUTCARファイルが正しい形式であるか、または存在するか確認してください。")