"""
VASP計算の出力ファイルを解析し、関連情報を提供するユーティリティモジュール。
このモジュールは、VASP (Vienna Ab initio Simulation Package) の計算結果ファイル
(INCAR, POSCAR, OUTCAR, EIGENVAL, DOSCAR, vasprun.xmlなど) を読み込み、
解析するための `tkVASP` クラスを提供します。
結晶構造情報、エネルギーデータ、バンド構造、状態密度、K点パスなどの情報を
抽出・整理し、Pythonオブジェクトとしてアクセス可能にします。
関連リンク: :doc:`tkvasp_usage`
"""
import os
import sys
from pathlib import Path
import shutil
import copy
import glob
import re
import csv
from pprint import pprint
import numpy as np
from numpy import sin, cos, tan, arccos, arcsin, arctan, sqrt, exp, log, pi
#import xml.etree.ElementTree as ET
#import xml.dom.minidom as md
import xmltodict
from pymatgen.io.vasp.outputs import Vasprun
from tklib.tkobject import tkObject
from tklib.tkutils import IsDir, IsFile, SplitFilePath, pconv, pint, pfloat, safe_getelement
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tkutils import SplitFilePath, GetList
from tklib.tkre import DelQuote, DelSpace
import tklib.tkre as tkre
from tklib.tkfile import tkFile
from tklib.tkcrystal.tkcrystalobject import tkCrystalObject
from tklib.tkcrystal.tkspacegroup import tkSpaceGroup
from tklib.tkprogvars import ProgramDir, DBDir, AtomDBDir, VASPDir, VASPPerlDir, VASPPythonDir
from tklib.tkfile import tkFile
from tklib.tkinifile import tkIniFile
from tklib.tksci.tksci import Reduce01, Round
from tklib.tksci.tksci import a0, NA, torad, todeg, Round, Factorize, IsFactor
from tklib.tksci.tksci import degcos, degsin, degtan, acos,asin, atan, degacos, degasin, degatan
from tklib.tksci.tkintegration import rieman_array_func
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
from tklib.tkcrystal.tkcrystal import tkCrystal
from tklib.tkvariousdata import tkVariousData
#from tklib.collabo.RH.parse_kpoints_opt import KPOINTS_OPT_Parser, write_kptopt_eig_csv, write_kptopt_dos_csv, write_high_symmetry_kpoints_opt
conv_dic = {
"alpha": 1
, "beta": 1
, "gamma": 1
, "delta": 1
, "epsilon": 1
, "zeta": 1
, "eta": 1
, "theta": 1
, "iota": 1
, "kappa": 1
, "lambda": 1
, "mu": 1
, "nu": 1
, "xi": 1
, "pi": 1
, "rho": 1
, "sigma": 1
, "tau": 1
, "upsilon": 1
, "phi": 1
, "chi": 1
, "psi": 1
, "omega": 1
}
[ドキュメント]
class tkVASP(tkCrystalObject):
"""
VASP計算の出力ファイルを解析し、関連情報を提供するクラス。
tkCrystalObjectを継承し、VASPのINCAR、POSCAR、POTCAR、KPOINTS、OUTCAR、EIGENVAL、
DOSCAR、vasprun.xmlなどの主要な出力ファイルからデータを読み込み、
結晶構造、エネルギー情報、バンド構造、状態密度などの解析機能を提供します。
"""
def __init__(self, **args):
"""
tkVASPクラスのコンストラクタ。
内部状態（結晶構造情報、INCAR情報など）を初期化します。
:param args: dict: キーワード引数として渡される初期設定。
"""
self.crystal = None
self.incarinf = None
self.kpointsinf = None
self.outcarinf = None
self.eigenvalinf = None
self.update(**args)
def __del__(self):
"""
tkVASPオブジェクトが破棄される際に呼び出されます。
現在は何も処理を行いません。
"""
pass
def __str__(self):
"""
オブジェクトの文字列表現を返します。
クラスのパスを文字列として返します。
:returns: str: クラスの完全パス。
"""
return self.ClassPath()
def __del__(self):
"""
tkVASPオブジェクトが破棄される際に呼び出されます。
現在は何も処理を行いません。
"""
pass
[ドキュメント]
def getdir(self, path):
"""
指定されたパスがディレクトリであればそのパスを、ファイルであればその親ディレクトリを返します。
ファイルパスからディレクトリ部分を抽出するために使用されます。
:param path: str: ファイルまたはディレクトリのパス。
:returns: str: 解決されたディレクトリのパス。
"""
if IsDir(path):
return path
dirname, basename, filebody, ext = SplitFilePath(path)
return dirname
[ドキュメント]
def get_CAR_dir(self, path, use_abspath = False):
"""
VASP計算の「CAR」ディレクトリのパスを取得します。
指定されたパスがディレクトリの場合はそのパスを、ファイルの場合はその親ディレクトリを返します。
`use_abspath` がTrueの場合、絶対パスに変換します。
:param path: str: ファイルまたはディレクトリのパス。
:param use_abspath: bool: 結果を絶対パスで返すかどうか。デフォルトはFalse。
:returns: str: VASP計算ディレクトリのパス。
"""
if use_abspath:
path = os.path.abspath(path)
if os.path.isdir(path):
return path
dirname = os.path.dirname(path)
return dirname
return os.path.join(CAR_dir, 'INCAR')
[ドキュメント]
def get_VASPPath(self, path, filename):
"""
VASP計算ディレクトリ内の特定のファイルの完全パスを取得します。
`get_CAR_dir` で取得したディレクトリとファイル名を結合します。
:param path: str: VASP計算ディレクトリ内にあるファイルまたはディレクトリのパス。
:param filename: str: 取得したいVASPファイルの名前（例: 'INCAR'）。
:returns: str: 指定されたVASPファイルの完全パス。
"""
CAR_dir = self.get_CAR_dir(path)
return os.path.join(CAR_dir, filename)
[ドキュメント]
def get_INCAR(self, CAR_dir):
"""
INCARファイルのパスを取得します。
`get_VASPPath` を使用してINCARファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: INCARファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'INCAR')
# return os.path.join(CAR_dir, 'INCAR')
[ドキュメント]
def get_POSCAR(self, CAR_dir):
"""
POSCARファイルのパスを取得します。
`get_VASPPath` を使用してPOSCARファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: POSCARファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'POSCAR')
# return os.path.join(CAR_dir, 'POSCAR')
[ドキュメント]
def get_POTCAR(self, CAR_dir):
"""
POTCARファイルのパスを取得します。
`get_VASPPath` を使用してPOTCARファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: POTCARファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'POTCAR')
# return os.path.join(CAR_dir, 'POTCAR')
[ドキュメント]
def get_KPOINTS(self, CAR_dir):
"""
KPOINTSファイルのパスを取得します。
`get_VASPPath` を使用してKPOINTSファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: KPOINTSファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'KPOINTS')
# return os.path.join(CAR_dir, 'KPOINTS')
[ドキュメント]
def get_KPOINTS_OPT(self, CAR_dir):
"""
KPOINTS_OPTファイルのパスを取得します。
`get_VASPPath` を使用してKPOINTS_OPTファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: KPOINTS_OPTファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'KPOINTS_OPT')
[ドキュメント]
def get_XML(self, CAR_dir):
"""
vasprun.xmlファイルのパスを取得します。
`get_VASPPath` を使用してvasprun.xmlファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: vasprun.xmlファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'vasprun.xml')
[ドキュメント]
def get_VASPRUN(self, CAR_dir):
"""
vasprun.xmlファイルのパスを取得します。
`get_VASPPath` を使用してvasprun.xmlファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: vasprun.xmlファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'vasprun.xml')
[ドキュメント]
def get_CONTCAR(self, CAR_dir):
"""
CONTCARファイルのパスを取得します。
`get_VASPPath` を使用してCONTCARファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: CONTCARファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'CONTCAR')
# return os.path.join(CAR_dir, 'CONTCAR')
[ドキュメント]
def get_OUTCAR(self, CAR_dir):
"""
OUTCARファイルのパスを取得します。
`get_VASPPath` を使用してOUTCARファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: OUTCARファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'OUTCAR')
# return os.path.join(CAR_dir, 'OUTCAR')
[ドキュメント]
def get_DOSCAR(self, CAR_dir):
"""
DOSCARファイルのパスを取得します。
`get_VASPPath` を使用してDOSCARファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: DOSCARファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'DOSCAR')
# return os.path.join(CAR_dir, 'DOSCAR')
[ドキュメント]
def get_EIGENVAL(self, CAR_dir):
"""
EIGENVALファイルのパスを取得します。
`get_VASPPath` を使用してEIGENVALファイルのパスを構築します。
:param CAR_dir: str: VASP計算ディレクトリのパス。
:returns: str: EIGENVALファイルの完全パス。
"""
return self.get_VASPPath(CAR_dir, 'EIGENVAL')
# return os.path.join(CAR_dir, 'EIGENVAL')
[ドキュメント]
def read_potcar_inf(self, potcarpath, print_level = 1):
"""
POTCARファイルから情報を読み込みます。
POTCARファイルの内容を解析し、各擬ポテンシャルに関する情報
（PP名、原子の種類、価電子数など）を辞書形式で抽出します。
:param potcarpath: str: POTCARファイルへのパス。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict or None: POTCAR情報を含む辞書、またはファイルが見つからない場合はNone。
"""
inf = {}
fp = tkFile(potcarpath, 'r')
if not fp or fp.fp is None:
return None
count = 0
pp_list = []
atom_list = []
Ne_list = []
while 1:
line = fp.ReadLine()
if not line:
break
# PAW_PBE O 08Apr2002
line = line.strip()
ret = tkre.Search(r'\S+\s+([A-Z][a-z]?)', line)
if ret:
pp = line.strip()
atom = ret[1]
line = fp.ReadLine()
Ne = pfloat(line)
inf["PP{}".format(count)] = pp
inf["Atom{}".format(count)] = ret[1]
inf["Ne{}".format(count)] = Ne
pp_list.append(pp)
atom_list.append(atom)
Ne_list.append(Ne)
fp.SkipTo("End of Dataset")
count += 1
inf["nPPTypes"] = count
inf["PPs"] = pp_list
inf["Atoms"] = atom_list
inf["Nes"] = Ne_list
self.potcarinf = inf
return inf
[ドキュメント]
def read_poscar_inf(self, poscarpath, is_print = None, print_level = 1):
"""
POSCARファイルから結晶構造情報を読み込みます。
POSCARファイルの内容を解析し、格子定数、格子ベクトル、原子の種類と数、
原子サイトの座標などを抽出します。`tkCrystal` オブジェクトを構築し、
内部状態として保持します。POTCARファイルも参照して原子タイプを特定する場合があります。
:param poscarpath: str: POSCARファイルへのパス。
:param is_print: bool or None: 結果を標準出力に表示するかどうか。Trueで表示、Falseで非表示。Noneの場合、`print_level`に従う。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict or None: POSCAR情報を含む辞書、またはファイルが見つからない場合はNone。
"""
if is_print is False:
print_level = 0
elif is_print is True:
print_level = 1
is_print = None
inf = {}
infp = tkFile(poscarpath, 'r')
if not infp or infp.fp is None:
return None
sample_name = infp.ReadLine().strip()
if print_level:
print("sample name: [{}]".format(sample_name))
a0 = pfloat(infp.ReadLine())
aij = []
aij.append(infp.ReadLine().split())
aij.append(infp.ReadLine().split())
aij.append(infp.ReadLine().split())
for i in range(3):
for j in range(3):
aij[i][j] = pfloat(aij[i][j])
if print_level:
print("")
print("Base lattice parameter: {} A".format(a0))
print("Lattice vectors:")
for i in range(3):
print(" {:16.12f} {:16.12f} {:16.12f}".format(*aij[i]))
line = infp.ReadLine()
# print("line=", line)
if re.search(r'^[\s\d\r\n]+$', line):
POTCAR_path = self.get_POTCAR(poscarpath)
# print("POTCAR_path=", POTCAR_path)
inf = self.read_potcar_inf(POTCAR_path)
# print("inf=", inf)
npp = inf['nPPTypes']
atomtypes = []
for i in range(npp):
atomtypes.append(inf.get(f'Atom{i}', None))
natoms = line.split()
else:
atomtypes = line.split()
natoms = infp.ReadLine().split()
for i in range(len(natoms)):
natoms[i] = pint(natoms[i])
if print_level:
print("")
print("Atom types and numbers:")
for i in range(len(atomtypes)):
print(" {:>4}: {:4d}".format(atomtypes[i], natoms[i]))
sitetypes = []
for ia in range(len(atomtypes)):
for n in range(natoms[ia]):
sitetypes.append(atomtypes[ia])
#Selective dynamics
line = infp.ReadLine().strip()
#Direct
if line.lower() != 'direct':
line = infp.ReadLine()
if print_level:
print("")
print("Atom positions:")
pos = []
count = 0
while 1:
prevpos = infp.Tell()
line = infp.ReadLine()
if not line:
break
pos0 = line.split()
if len(pos0) < 3:
break
try:
pos0[0] = float(pos0[0])
except:
currpos = infp.Seek(prevpos)
break
pos.append([pfloat(pos0[0]), pfloat(pos0[1]), pfloat(pos0[2])])
if print_level:
print(" {:4}: ({:16.12f}, {:16.12f}, {:16.12f})"
.format(sitetypes[count], pos[count][0], pos[count][1], pos[count][2]))
count += 1
if print_level:
print("")
print("Build Crystal object:")
cry = tkCrystal()
cry.SetSampleName(sample_name)
cry.SetCrystalName(sample_name)
aij = a0 * np.array(aij)
# print("aij_original=\n", aij)
cry.SetLatticeVectors(aij)
# print("aij_cry=\n", cry.LatticeVectors())
latt = cry.LatticeParameters()
sites = []
for i in range(len(pos)):
cry.AddAtomSite(name = sitetypes[i], pos = pos[i])
sites.append({ "atom": sitetypes[i], "x": pos[i][0], "y": pos[i][1], "z": pos[i][2]})
cry.ExpandCoordinates()
# cry.PrintInf()
inf["sample_name"] = sample_name
inf["crystal"] = cry
inf["a0"] = a0
inf["aij"] = np.array(a0) * aij
inf["nAtomTypes"] = len(atomtypes)
inf["AtomTypes"] = atomtypes
inf["nAtomSites"] = len(sites)
inf["AtomSites"] = sites
inf["LatticeParameters"] = [float(v) for v in latt]
inf["Vcell"] = cry.Volume()
self.crystal = cry
self.poscarinf = inf
return inf
[ドキュメント]
def read_poscar(self, poscarpath, is_print = None, print_level = 1):
"""
POSCARファイルから `tkCrystal` オブジェクトを読み込みます。
`read_poscar_inf` メソッドを呼び出し、その結果から `tkCrystal` オブジェクトを返します。
:param poscarpath: str: POSCARファイルへのパス。
:param is_print: bool or None: 結果を標準出力に表示するかどうか。Trueで表示、Falseで非表示。Noneの場合、`print_level`に従う。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: tkCrystal or None: 構築された `tkCrystal` オブジェクト、またはファイルが見つからない場合はNone。
"""
inf = self.read_poscar_inf(poscarpath, is_print = is_print, print_level = print_level)
return inf["crystal"]
[ドキュメント]
def read_incar_inf(self, incarpath, print_level = 1):
"""
INCARファイルから情報を読み込みます。
INCARファイルの内容を解析し、VASP計算のパラメータ（SYSTEM、各種設定値など）を
辞書形式で抽出します。コメント行や空白は無視されます。
:param incarpath: str: INCARファイルへのパス。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict or None: INCAR情報を含む辞書、またはファイルが見つからない場合はNone。
"""
fp = tkFile(incarpath, 'r')
if not fp or fp.fp is None:
return None
inf = {'SYSTEM': ''}
while 1:
line = fp.ReadLine()
if not line:
break
line = tkre.Sub(r'[!#].*$', '', line)
# print("line=", line)
aa = line.split(';')
for i in range(len(aa)):
ret = fp.Search(r'(\S+)\s*[=:]\s*(\S.*)\s*', aa[i])
if ret is None:
continue
# print("key=", ret)
inf[ret[1]] = pconv(ret[2])
self.incarinf = inf
return inf
[ドキュメント]
def read_oszicar_inf(self, oszicarpath):
"""
OSZICARファイルから情報を読み込みます。
OSZICARファイルの内容を解析し、各イオンスステップにおける電子ステップ、
全エネルギー、自由エネルギーなどの収束情報を抽出します。
:param oszicarpath: str: OSZICARファイルへのパス。
:returns: dict or None: OSZICAR情報を含む辞書、またはファイルが見つからない場合はNone。
"""
fp = tkFile(oszicarpath, 'r')
if not fp or fp.fp is None:
return None
inf = {}
ion_steps = []
all_steps = []
ion_step = 1
while 1:
line = fp.ReadLine()
if not line or line == '':
break
# print("l: ", line, end = '')
if 'DAV:' in line:
aa = line.split()
all_steps.append({ 'ion_step': ion_step, 'electron_step': pint(aa[1]), 'E': pfloat(aa[2]), 'dE': pfloat(aa[3])})
# print("E=", aa[2])
elif 'F=' in line:
aa = line.split()
ion_steps.append({ 'ion_step': aa[0], 'F': pfloat(aa[2]), 'E0': pfloat(aa[4]), 'dE': pfloat(aa[6])})
# print("F=", aa[2], " E=", aa[4])
ion_step += 1
elif 'ncg' in line:
pass
else:
pass
fp.Close()
inf['ion_steps'] = ion_steps
inf['all_steps'] = all_steps
self.oszicarinf = inf
return inf
[ドキュメント]
def read_outcar_inf(self, outcarpath, print_level = 1):
"""
OUTCARファイルから計算結果情報を読み込みます。
OUTCARファイルの内容を詳細に解析し、VASPの実行環境、INCARパラメータ、
フェルミエネルギー、平均コアポテンシャル、総エネルギー履歴、最終的な電荷、磁化、
誘電率、圧電係数、有効Born電荷、収束情報などを辞書形式で抽出します。
:param outcarpath: str: OUTCARファイルへのパス。
:returns: dict or None: OUTCAR情報を含む辞書、またはファイルが見つからない場合はNone。
"""
fp = tkFile(outcarpath, 'r')
if not fp or fp.fp is None: return None
inf = { "binary": "", "compiler": "", "mpi": "", "ncores_k": "", "ncores_band": ""}
line = fp.ReadLine()
if 'vasp' in line:
inf["binary"] = line.strip()
line = fp.SkipTo(r'^\s*executed', max_lines = 20)
if line is not None:
inf["compiler"] = line.replace("executed on", "").strip()
line = fp.SkipTo(r'^\s*running', max_lines = 20)
if line:
inf["mpi"] = line.replace("running", "").strip()
line = fp.SkipTo(r'^\s*distrk:', max_lines = 20)
inf["ncores_k"] = line.replace("distrk:", "").strip()
line = fp.SkipTo(r'^\s*distr:', max_lines = 20)
inf["ncores_band"] = line.replace("distr:", "").strip()
fp.Rewind()
line = fp.skip_to(r'^\s*INCAR:')
while 1:
line = fp.ReadLine()
if not line: break
# if '----------' in line:
if 'POTCAR:' in line: break
# print("line=", line)
ret = tkre.Search(r'(\S+)\s*=\s*(\S+)', line)
if ret is not None and len(ret) >= 3:
inf[ret[1]] = pconv(ret[2])
# print(f"{ret[1]}={ret[2]}")
line = fp.SkipTo("Dimension of arrays:")
while True:
line = fp.readline()
if not line: break
line = line.strip()
if line == "": break
pattern = r"(\w+)\s*=\s*(\d+)"
matches = re.findall(pattern, line)
for key, value in matches:
inf[key] = pconv(value)
inf["SYSTEM"] = ""
ret = fp.search_to(r'SYSTEM\s*=\s*(.*)\s*$')
if ret and len(ret) > 1: inf["SYSTEM"] = ret[1]
line = fp.SkipTo("Startparameter for this run:")
while 1:
line = fp.ReadLine()
if not line: break
if '----------' in line: break
aa = line.split(';')
for i in range(len(aa)):
ret = tkre.Search(r'(\S+)\s*=\s*(\S+)', aa[i])
if ret is None: continue
if len(ret) >= 3:
inf[ret[1]] = pconv(ret[2])
m = fp.SearchTo(r"E-fermi :\s*([+\-\d\.eDdD]+)", defval = "0.0")
EF = pconv(m[1])
inf["EF"] = EF
inf["E-fermi"] = EF
fp.rewind()
# Read average core potnetials
Vcore = []
# average (electrostatic) potential at core
line = fp.skip_to("potential at core")
# the test charge radii are 1.2793 0.7215
fp.ReadLine()
# (the norm of the test charge is 1.0000)
fp.ReadLine()
# 1 -52.1673 2 -52.1683 3 -52.1600 4 -52.1667 5 -70.7270
# 6 -70.7333 7 -70.7334 8 -70.7355 9 -70.7342 10 -70.7302
icount = 0
while 1:
line = fp.ReadLine()
if not line:
break
aa = line.split()
if len(aa) <= 1:
break
for i in range(0, len(aa), 2):
Vcore.append(pfloat(aa[i+1]))
inf["Vcore"] = Vcore
fp.rewind()
# read ion step variation of converged total energy
count = 0
energy = []
energy_without_entropy = []
energy_zero_sigma = []
while 1:
line = fp.SkipTo("\\s*free energy TOTEN =");
if not line:
break
pos = fp.Tell()
s = tkre.Search(r'free energy TOTEN =\s*(\S+)', line)
# if pos > 0:
# fp.Seek(pos, 0)
key = 'FreeEnergy{}'.format(count)
energy.append(pfloat(s[1]))
line = fp.SkipTo("\\s*energy without entrop");
if not line:
break
s = tkre.Search(r'=\s*(\S+)\s.*=\s*(\S+)', line)
energy_without_entropy.append(pfloat(s[1]))
energy_zero_sigma.append(pfloat(s[2]))
if pos > 0:
fp.Seek(pos, 0)
# print("E=", energy[count])
count += 1
inf['FreeEnergy'] = energy
inf['TOTEN'] = energy[len(energy) - 1]
inf['TOTEN_without_entropy'] = energy_without_entropy[len(energy_without_entropy) - 1]
inf['TOTEN_zero_sigma'] = energy_zero_sigma[len(energy_zero_sigma) - 1]
# read final charges
line = fp.SkipTo("\\s*total charge", 0)
line = fp.SkipTo("# of ion")
line = fp.ReadLine()
def split_spdf(_a):
na = len(_a)
if na >= 1:
iIon = pint(_a[0])
for i in range(1, na):
_a[i] = pfloat(_a[i])
if na == 6:
s, p, d, f, tot = _a[1], _a[2], _a[3], _a[4], _a[5]
elif na == 5:
s, p, d, f, tot = _a[1], _a[2], _a[3], 0.0, _a[4]
elif na == 4:
s, p, d, f, tot = _a[1], _a[2], 0.0, 0.0, _a[3]
elif na == 3:
s, p, d, f, tot = _a[1], 0.0, 0.0, 0.0, _a[2]
return { 's'.format(iIon): s,
'p'.format(iIon): p,
'd'.format(iIon): d,
'f'.format(iIon): f,
'tot'.format(iIon): tot
}
count = 0
chargeinf = []
while 1:
line = fp.ReadLine()
# print("l", line.strip())
if not line: break
line = line.strip()
if tkre.Match(r'^-----', line): break
_a = line.split()
chargeinf.append(split_spdf(_a))
count += 1
# print("line=", line)
inf['nIon'] = ""
inf['FinalCharges'] = []
inf['FinalTotalCharges'] = []
inf['FinalMagnetization'] = []
inf['FinalTotalMagnetizations'] = []
if line:
inf['FinalCharges'] = chargeinf
inf['nIon'] = count
line = fp.readline()
# print("line=", line)
_a = line.split()
inf['FinalTotalCharges'] = split_spdf(_a)
# read final charges
line = fp.SkipTo("\\s*magnetization", 0)
if line:
line = fp.SkipTo("# of ion")
line = fp.ReadLine()
count = 0
spininf = []
while 1:
line = fp.ReadLine()
if not line: break
line = line.strip()
if tkre.Match(r'^-----', line): break
_a = line.split()
spininf.append(split_spdf(_a))
count += 1
line = fp.readline()
_a = line.split()
inf['FinalTotalMagnetizations'] = split_spdf(_a)
inf['FinalMagnetization'] = spininf
pos = fp.tell()
# Check EF after charge lines
# BZINTS: Fermi energy: 2.623512;********** electrons
fp.Rewind()
while 1:
m = fp.SearchTo(r"BZINTS: Fermi energy:\s*([+\-\d\.eDdD]+)", defval = None)
if m is None:
break
else:
EF = pconv(m[1])
inf["EF"] = EF
inf["EF(BZINTS)"] = EF
# if 'T' in sainf['LHFCALC']:
if 'T' in inf.get('LHFCALC', 'F'):
inf['IsHF'] = 1
else:
inf['IsHF'] = 0
if inf.get('METAGGA', '') != '':
inf['IsMETAGGA'] = 1
else:
inf['IsMETAGGA'] = 0
# fp.seek(pos, 0)
fp.rewind()
eps_static_localeff = None
line = fp.SkipTo(r"MACROSCOPIC STATIC DIELECTRIC TENSOR \(including local field effects")
if line:
eps_static_localeff = []
fp.readline()
_a = fp.readline().split()
eps_static_localeff.append([float(v) for v in _a])
_a = fp.readline().split()
eps_static_localeff.append([float(v) for v in _a])
_a = fp.readline().split()
eps_static_localeff.append([float(v) for v in _a])
inf["eps_static_localeff"] = eps_static_localeff
piezo_static_localeff = None
line = fp.SkipTo(r"PIEZOELECTRIC TENSOR \(including local field effects\) for field in x, y, z\s*\(C")
if line:
piezo_static_localeff = []
fp.readline()
fp.readline()
_a = fp.readline().split()
piezo_static_localeff.append([float(v) for v in _a[1:]])
_a = fp.readline().split()
piezo_static_localeff.append([float(v) for v in _a[1:]])
_a = fp.readline().split()
piezo_static_localeff.append([float(v) for v in _a[1:]])
inf["piezo_static_localeff"] = piezo_static_localeff
born_charges = None
line = fp.SkipTo(r"BORN EFFECTIVE CHARGES \(including local field effects")
if line:
born_charges = []
fp.readline()
while True:
line = fp.readline()
if not line or line.strip() == '': break
_a = line.split()
iion = int(_a[1])
born_charges.append([])
_a = fp.readline().split()
born_charges[iion-1].append([float(v) for v in _a])
_a = fp.readline().split()
born_charges[iion-1].append([float(v) for v in _a])
_a = fp.readline().split()
born_charges[iion-1].append([float(v) for v in _a])
inf["born_charges"] = born_charges
eps_static_ionic = None
line = fp.SkipTo(r"MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC CONTRIBUTION")
if line:
eps_static_ionic = []
fp.readline()
_a = fp.readline().split()
eps_static_ionic.append([float(v) for v in _a])
_a = fp.readline().split()
eps_static_ionic.append([float(v) for v in _a])
_a = fp.readline().split()
eps_static_ionic.append([float(v) for v in _a])
inf["eps_static_ionic"] = eps_static_ionic
fp.rewind()
convergence_list = []
while True:
line = fp.skip_to("aborting loop because EDIFF is reached")
if not line: break
convergence_list.append(line.strip())
inf["EDIFF convergence"] = convergence_list
fp.rewind()
inf["EDIFFG convergence"] = ""
line = fp.skip_to("reached required accuracy")
if line:
inf["EDIFFG convergence"] = line.strip()
# fp.seek(pos, 0)
fp.rewind()
for key in ["Total CPU time", "User time", "System time", "Elapsed time"]:
m = fp.SearchTo(rf"^\s*{key}.*:\s*([0-9\.]+)", defval = None)
if m:
inf[key] = pfloat(m[1])
fp.Close()
self.outcarinf = inf
return inf
[ドキュメント]
def read_n_crystalstructures(self, outcarpath):
"""
OUTCARファイルから収束した結晶構造の数を読み込みます。
OUTCARファイル内で「aborting loop because EDIFF is reached」という文字列が
出現する回数を数えることで、構造最適化の各ステップで収束した構造の数を特定します。
:param outcarpath: str: OUTCARファイルへのパス。
:returns: int: 収束した結晶構造の数。ファイルが見つからない場合は0。
"""
fp = tkFile(outcarpath, 'r')
n = 0
if not fp:
return 0
while 1:
line = fp.SkipTo('aborting loop because EDIFF is reached')
if not line:
break
n += 1
return n
[ドキュメント]
def read_next_crystalstructure(self, fp, template_cry, IsReduce01 = False):
"""
OUTCARファイルから次の結晶構造と関連情報を読み込みます。
`read_n_crystalstructures` で取得される各構造最適化ステップの情報を順次読み取ります。
格子ベクトル、原子座標、原子にかかる力、総エネルギーなどを抽出し、`tkCrystal` オブジェクトを更新します。
:param fp: tklib.tkfile.tkFile: OUTCARファイルのファイルポインタ。
:param template_cry: tkCrystal: テンプレートとなる `tkCrystal` オブジェクト（原子種情報などを含む）。
:param IsReduce01: bool: 原子座標を0-1の範囲に正規化するかどうか。デフォルトはFalse。
:returns: tuple: (dict, tkCrystal) または (None, None): 抽出された情報辞書と更新された `tkCrystal` オブジェクト、または次の構造が見つからない場合は(None, None)。
"""
cry = copy.deepcopy(template_cry)
inf = {}
line = fp.SkipTo('aborting loop because EDIFF is reached')
if not line:
return None, None
line = fp.SkipTo('VOLUME and BASIS-vectors')
for i in range(4):
fp.ReadLine()
a1 = fp.ReadLine().split()
a2 = fp.ReadLine().split()
a3 = fp.ReadLine().split()
a1[0], a1[1], a1[2] = pfloat(a1[0]), pfloat(a1[1]), pfloat(a1[2])
a2[0], a2[1], a2[2] = pfloat(a2[0]), pfloat(a2[1]), pfloat(a2[2])
a3[0], a3[1], a3[2] = pfloat(a3[0]), pfloat(a3[1]), pfloat(a3[2])
cry.SetLatticeVectors([a1[:3], a2[:3], a3[:3]])
latt = cry.LatticeParameters()
line = fp.SkipTo('POSITION')
# skip ' ------...'
line = fp.ReadLine()
count = 0
while 1:
line = fp.ReadLine()
if not line or '-----' in line:
break
x, y, z, fx, fy, fz, = line.split()
# x = pfloat(x) / latt[0]
# y = pfloat(y) / latt[1]
# z = pfloat(z) / latt[2]
# Coordinate
frx, fry, frz = cry.CartesianToFractional(pfloat(x), pfloat(y), pfloat(z))
if IsReduce01:
frx, fry, frz = Reduce01(frx), Reduce01(fry), Reduce01(frz)
# print("")
# print("***xyz=({},{},{})".format(x,y,z))
# print(" frxyz=({},{},{})".format(frx,fry,frz))
# x2,y2,z2 = cry.FractionalToCartesian(frx, fry, frz)
# print(" 2xyz=({},{},{})".format(x2,y2,z2))
# exit()
# Force
fx, fy, fz = cry.CartesianToFractional(pfloat(fx), pfloat(fy), pfloat(fz))
cry.SetAtomSite(count, pos = [frx, fry, frz],
force = [pfloat(fx), pfloat(fy), pfloat(fz)])
count += 1
cry.ExpandCoordinates()
line = fp.SkipTo("\\s*free energy TOTEN =");
ret = fp.Search(r'free energy TOTEN =\s*(\S+)', line)
inf["TOTEN"] = pfloat(ret[1])
pos = fp.Tell()
# line = fp.SkipTo('MD-specific')
line = fp.SkipTo('EKIN')
if line:
ret = fp.Search(r'EKIN\s*=\s*([+\-\d\.eEdD]+)', line, defval = ['', ''])
inf["EKIN"] = pfloat(ret[1])
line = fp.SkipTo('temperature')
ret = fp.Search(r'EKIN_LAT\s*=\s*([+\-\d\.eEdD]+)', line, defval = ['', ''])
inf["EKIN_LAT"] = pfloat(ret[1])
ret = fp.Search(r'temperature\s+([+\-\d\.eEdD]+)', line, defval = ['', ''])
inf["temperature"] = pfloat(ret[1])
line = fp.SkipTo('ETOTAL')
ret = fp.Search(r'ETOTAL\s*=\s*([+\-\d\.eEdD]+)', line, defval = ['', ''])
inf["ETOTAL"] = pfloat(ret[1])
else:
fp.Seek(pos)
return inf, cry
[ドキュメント]
def convert_kname(self, name):
"""
K点パスの記号名をLaTeX形式に変換します。
K点パスの記号名（例: "Gamma"）を認識し、対応するLaTeX形式（例: "$\\Gamma$"）に変換します。
ギリシャ文字などが対象です。
:param name: str: 変換するK点パスの記号名。
:returns: str: 変換された記号名、または変換対象外の場合は元の名前。
"""
if len(name) <= 1:
return name
try:
a = conv_dic[name.lower()]
name = name.capitalize()
name = '$\\{}$'.format(name)
except:
pass
return name
[ドキュメント]
def read_kpoints_linemode(self, fp, inf, cry = None):
"""
KPOINTSファイル（ラインモード）からK点パス情報を読み込みます。
バンド計算などで使用されるラインモードのKPOINTSファイルを解析し、
K点の座標、重み、K点間の距離、パス上の累積距離などを抽出します。
:param fp: tklib.tkfile.tkFile: KPOINTSファイルのファイルポインタ。
:param inf: dict: 情報を格納する辞書。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。K点距離計算に必要。
:returns: dict: 更新されたK点パス情報を含む辞書。
"""
IsEOF = 0
kp = []
kp_scf = []
kp_all = []
kpdic = []
ktot = 0.0
while 1:
while 1:
line1 = fp.ReadLine()
if not line1:
IsEOF = 1
break
line1 = line1.strip()
if line1 == '':
next
else:
break
if IsEOF:
break
while 1:
line2 = fp.ReadLine()
if not line2:
IsEOF = 1
break
line2 = line2.strip()
if line2 == '':
next
else:
break
if IsEOF:
break
# print("line1: ", line1)
# print("line2: ", line2)
a1 = line1.split()
n1 = len(a1)
kx0, ky0, kz0 = pfloat(a1[0]), pfloat(a1[1]), pfloat(a1[2])
a2 = line2.split()
n2 = len(a2)
kx1, ky1, kz1 = pfloat(a2[0]), pfloat(a2[1]), pfloat(a2[2])
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx1, ky1, kz1])
ktot += dk
match = tkre.Search(r'[!#]\s*(\S+)', line1)
if match:
name0 = match[1]
else:
name0 = ''
match = tkre.Search(r'[!#]\s*(\S+)', line2)
if match:
name1 = match[1]
else:
name1 = ''
name0c = self.convert_kname(name0)
name1c = self.convert_kname(name1)
w = 0.0
kinf = {"kx0": kx0, "ky0": ky0, "kz0": kz0, "kname0": name0, "kname0_conv": name0c,
"kx1": kx1, "ky1": ky1, "kz1": kz1, "kname1": name1, "kname1_conv": name1c,
"kw": w, "dk": dk, "ktot": ktot}
kp.append(kinf)
kp_all.append(kinf.copy())
kpdic.append(kinf.copy())
inf['nk'] = len(kp)
inf['kpoints'] = kp
inf['kpoints_scf'] = kp_scf
inf['kpoints_all'] = kp_all
inf['kpointsdic'] = kpdic
return inf
[ドキュメント]
def read_kpoints_pointmode(self, fp, inf, cry = None):
"""
KPOINTSファイル（ポイントモード）からK点情報を読み込みます。
自己無撞着計算などで使用されるポイントモードのKPOINTSファイルを解析し、
K点の座標、重み、K点間の距離、パス上の累積距離などを抽出します。
特定のK点（重みが0のK点）は別途抽出されます。
:param fp: tklib.tkfile.tkFile: KPOINTSファイルのファイルポインタ。
:param inf: dict: 情報を格納する辞書。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。K点距離計算に必要。
:returns: dict: 更新されたK点情報を含む辞書。
"""
kp = []
kp_scf = []
kp_all = []
while 1:
line1 = fp.ReadLine()
if not line1:
break
# print("line1: ", line1)
a1 = line1.split()
n1 = len(a1)
if n1 < 4:
continue
kx1, ky1, kz1, w = pfloat(a1[0]), pfloat(a1[1]), pfloat(a1[2]), pfloat(a1[3])
match = tkre.Search(r'[!#]\s*(\S+)', line1)
if match:
name1 = match[1]
else:
name1 = ''
name1c = self.convert_kname(name1)
dk = None
ktot = None
kinf = {"kx0": kx1, "ky0": ky1, "kz0": kz1, "kname0": name1, "kname0_conv": name1c,
"kx1": kx1, "ky1": ky1, "kz1": kz1, "kname1": name1, "kname1_conv": name1c,
"kw": w, "dk": dk, "ktot": ktot}
kp_all.append(kinf)
if w == 0.0:
kp.append(kinf.copy())
else:
kp_scf.append(kinf.copy())
kx0, ky0, kz0 = kx1, ky1, kz1
nk = len(kp)
kpdic = []
ktot = 0.0
dkdic = 0.0
dkvec0 = None
for ik in range(nk):
ki = kp[ik]
kx1, ky1, kz1, w = ki["kx0"], ki["ky0"], ki["kz0"], ki["kw"]
if ik == 0:
kx0, ky0, kz0 = kx1, ky1, kz1
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx1, ky1, kz1])
ktot += dk
dkdic += dk
ki["dk"] = dk
ki["ktot"] = ktot
# add to kpdic
if ik == 0:
kpdic.append(ki.copy())
elif ik == nk - 1:
kpdic.append(ki.copy())
else:
kim = kp[ik-1]
kip = kp[ik+1]
kx0, ky0, kz0, w = kim["kx0"], kim["ky0"], kim["kz0"], kim["kw"]
kx2, ky2, kz2, w = kip["kx0"], kip["ky0"], kip["kz0"], kip["kw"]
r2 = (kx1 - kx0) * (kx1 - kx0) + (ky1 - ky0) * (ky1 - ky0) + (kz1 - kz0) * (kz1 - kz0)
innp = (kx1 - kx0) * (kx2 - kx1) + (ky1 - ky0) * (ky2 - ky1) + (kz1 - kz0) * (kz2 - kz1)
if r2 > 1.0e-6 and (name1 != '' or abs(innp) < 1.0e-6):
kidic = ki.copy()
kidic["dk"] = dkdic
kpdic.append(kidic)
dkdic = 0.0
kx0, ky0, kz0 = kx1, ky1, kz1
ktot = 0.0
for ik in range(len(kp_scf)):
ki = kp_scf[ik]
kx1, ky1, kz1, w = ki["kx0"], ki["ky0"], ki["kz0"], ki["kw"]
if ik == 0:
kx0, ky0, kz0 = kx1, ky1, kz1
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx1, ky1, kz1])
ktot += dk
ki["dk"] = dk
ki["ktot"] = ktot
kx0, ky0, kz0 = kx1, ky1, kz1
ktot = 0.0
for ik in range(len(kp_all)):
ki = kp_all[ik]
kx1, ky1, kz1, w = ki["kx0"], ki["ky0"], ki["kz0"], ki["kw"]
if ik == 0:
kx0, ky0, kz0 = kx1, ky1, kz1
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx1, ky1, kz1])
ktot += dk
ki["dk"] = dk
ki["ktot"] = ktot
kx0, ky0, kz0 = kx1, ky1, kz1
inf['nk'] = len(kp)
inf['kpoints'] = kp
inf['kpoints_scf'] = kp_scf
inf['kpoints_all'] = kp_all
inf['kpointsdic'] = kpdic
return inf
[ドキュメント]
def read_kpoints(self, KPOINTS_path, cry = None):
"""
KPOINTSファイルからK点情報を読み込みます。
KPOINTSファイルのモード（ラインモードかポイントモードか）を自動的に判別し、
それぞれに対応するパーサー (`read_kpoints_linemode` または `read_kpoints_pointmode`) を呼び出します。
:param KPOINTS_path: str: KPOINTSファイルへのパス。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。K点距離計算に必要。
:returns: dict or None: K点情報を含む辞書、またはファイルが見つからない場合はNone。
"""
if not os.path.isfile(KPOINTS_path):
return None
if cry is None:
cry = self.crystal
inf = {}
fp = tkFile(KPOINTS_path, 'r')
if fp is None:
return None
inf['title'] = fp.ReadLine().strip()
inf['nk_target'] = pint(fp.ReadLine())
mode = fp.ReadLine().strip().upper()
if mode[0] == 'L':
inf['mode'] = mode
inf['coord_mode'] = fp.ReadLine().strip()
inf = self.read_kpoints_linemode(fp, inf, cry)
else:
inf['mode'] = mode
inf['coord_mode'] = mode
inf = self.read_kpoints_pointmode(fp, inf, cry)
fp.Close()
self.kpointsinf = inf
return inf
[ドキュメント]
def read_eigenval_opt(self, EIGENVAL_path, normalize_E = True, KPOINTS_inf = None,
ISPIN = None, IsHF = None, IsMETAGGA = None, EF = None, NELECT = None, cry = None,
print_level = 1):
"""
`vasprun.xml` と `KPOINTS_OPT` からバンド構造情報を読み込みます。
`pymatgen` の `Vasprun` クラスを使用せず、`KPOINTS_OPT_Parser` を用いて、
`vasprun.xml` と `KPOINTS_OPT` からEIGENVAL情報（バンドエネルギーと占有数）を抽出します。
これは特にバンドパスが複雑な場合や特定のK点に重みがない場合のバンド構造解析に適しています。
:param EIGENVAL_path: str: EIGENVALファイルへのパス（KPOINTS_OPTとXMLのパスを特定するために使用）。
:param normalize_E: bool: エネルギーをフェルミ準位で正規化するかどうか。デフォルトはTrue。
:param KPOINTS_inf: dict or None: KPOINTS情報を含む辞書。指定されない場合、内部で読み込まれる。
:param ISPIN: int or None: スピン偏極計算 (`ISPIN=2`) かどうか。Noneの場合、OUTCAR情報から取得。
:param IsHF: bool or None: Hartree-Fock計算かどうか。Noneの場合、OUTCAR情報から取得。
:param IsMETAGGA: bool or None: meta-GGA計算かどうか。Noneの場合、OUTCAR情報から取得。
:param EF: float or None: フェルミエネルギー。Noneの場合、OUTCAR情報から取得。
:param NELECT: int or None: 総電子数。Noneの場合、OUTCAR情報から取得。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict: バンド構造情報を含む辞書。
"""
if self.outcarinf is None: return None
if cry is None:
cry = self.crystal
if NELECT is None:
NELECT = self.outcarinf["NELECT"]
if ISPIN is None:
ISPIN = self.outcarinf["ISPIN"]
if IsHF is None:
IsHF = self.outcarinf["IsHF"]
if IsMETAGGA is None:
IsMETAGGA = self.outcarinf["IsMETAGGA"]
if EF is None:
EF = self.outcarinf["EF"]
KPOINTS_OPT_path = self.get_KPOINTS_OPT(EIGENVAL_path)
XML_path = self.get_XML(EIGENVAL_path)
kpoints_inf = self.read_kpoints(KPOINTS_OPT_path, cry = cry)
# print("kpoints_inf=", kpoints_inf)
if kpoints_inf is None:
return {}
Nline = 20
parser = KPOINTS_OPT_Parser(XML_path, Nline)
EF = parser.get_fermi_level()
# parser.reciprocal_lattice = self.vasprun.final_structure.lattice.reciprocal_lattice.matrix
# parser.bands_data = self.get_eigenvalues()
# parser.kpath = self.get_kpath()
# bands_data = self.get_eigenvalues()
# bands_data.append({
# 'spin': spin,
# 'kpoint': kpoint,
# 'eigenvalue': float(r.text.strip()) - fermi_level
# })
kpoints = parser.kpoints
weights = parser.weights_list
kpath = parser.kpath
inf = {}
bands_data = parser.get_eigenvalues()
eigenvalues_dict, occupancies_dict = parser.eigenvalues_per_kpoint()
inf["n"] = None
inf["nk"] = len(kpoints)
inf["nLevels"] = max(len(eigenvalues) for eigenvalues in eigenvalues_dict.values())
# for i, k in enumerate(kpoints):
# print("k=", i, k)
# for i, eig in enumerate(eigenvalues_dict.values()):
# print("eig=", i, len(eig))
# exit()
# print("nk=", inf['nk'])
# print("nLevels=", inf['nLevels'])
ktot_prev = 0.0
Elist = [[]] * inf["nk"]
# Elist_scf = []
# Elist_all = []
kx0 = None
ky0 = None
kz0 = None
for key, E in eigenvalues_dict.items():
spin, ik = key
if '1' in spin:
spin = 0
else:
spin = 1
occ = occupancies_dict[key]
ktot = kpath[ik]
wk = weights[ik]
kv = kpoints[ik]
dk = ktot - ktot_prev
# print(f"{spin} {ik}: {ktot} {wk} {kv} {occ}: {E}")
if len(Elist[ik]) > 2:
_kv0, _kv1, _kv2, _wk, _dk, _ktot, _Eup, _occup, _Edn, _occdn = Elist
else:
_Eup = None
_occup = None
_Edn = None
_occdn = None
if spin == 0:
Elist[ik] = [kv[0], kv[1], kv[2], wk, dk, ktot, E, occ, _Edn, _occdn]
else:
Elist[ik] = [kv[0], kv[1], kv[2], wk, dk, ktot, _Eup, _occup, E, occ]
ktot_prev = ktot
# print()
# c = 0
# for E in Elist:
# print("E=", E)
# if len(E) < 2:
# c += 1
# print("nk=", inf["nk"])
# print("c=", c)
inf['EList'] = Elist
inf['EList_scf'] = Elist
inf['EList_all'] = Elist
self.eigenvalinf = inf
return inf
[ドキュメント]
def read_eigenval(self, EIGENVAL_path, normalize_E = True, KPOINTS_inf = None, ISPIN = None,
IsHF = None, IsMETAGGA = None, EF = None, NELECT = None, cry = None,
print_level = 1):
"""
EIGENVALファイルからバンド構造情報を読み込みます。
EIGENVALファイルの内容を解析し、各K点におけるバンドエネルギーと占有数を抽出します。
スピン偏極計算の場合、スピンアップとスピンダウンの情報を区別して読み込みます。
K点パスが定義されている場合（ラインモード）、K点間の距離も計算します。
:param EIGENVAL_path: str: EIGENVALファイルへのパス。
:param normalize_E: bool: エネルギーをフェルミ準位で正規化するかどうか。デフォルトはTrue。
:param KPOINTS_inf: dict or None: KPOINTS情報を含む辞書。指定されない場合、内部で読み込まれる。
:param ISPIN: int or None: スピン偏極計算 (`ISPIN=2`) かどうか。Noneの場合、OUTCAR情報から取得。
:param IsHF: bool or None: Hartree-Fock計算かどうか。Noneの場合、OUTCAR情報から取得。
:param IsMETAGGA: bool or None: meta-GGA計算かどうか。Noneの場合、OUTCAR情報から取得。
:param EF: float or None: フェルミエネルギー。Noneの場合、OUTCAR情報から取得。
:param NELECT: int or None: 総電子数。Noneの場合、OUTCAR情報から取得。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict or None: バンド構造情報を含む辞書、またはファイルが見つからない場合はNone。
"""
if self.outcarinf is None: return None
if not os.path.exists(EIGENVAL_path): return None
INCAR_path = self.get_INCAR(EIGENVAL_path)
incar = self.read_incar_inf(INCAR_path)
if incar:
LKPOINTS_OPT = incar.get('LKPOINTS_OPT', None)
if LKPOINTS_OPT is None or 'false' in LKPOINTS_OPT.lower():
LKPOINTS_OPT = False
else:
LKPOINTS_OPT = True
else:
LKPOINTS_OPT = True
KPOINTS_OPT_path = self.get_KPOINTS_OPT(EIGENVAL_path)
if LKPOINTS_OPT and os.path.isfile(KPOINTS_OPT_path):
return self.read_eigenval_opt(EIGENVAL_path, normalize_E, KPOINTS_inf, ISPIN, IsHF, IsMETAGGA, EF, NELECT, cry)
if cry is None:
cry = self.crystal
if NELECT is None:
NELECT = self.outcarinf["NELECT"]
if ISPIN is None:
ISPIN = self.outcarinf["ISPIN"]
if IsHF is None:
IsHF = self.outcarinf["IsHF"]
if EF is None:
EF = self.outcarinf["EF"]
try:
kptinf = KPOINTS_inf['kpoints_all']
except:
kptinf = None
# print("kpinf=", kpinf)
inf = {}
fp = tkFile(EIGENVAL_path, 'r')
if fp.fp is None: return None
# 4 4 1 1
fp.ReadLine()
# 0.1236685E+02 0.3291500E-09 0.3291500E-09 0.5272308E-09 0.5000000E-15
fp.ReadLine()
# 1.000000000000000E-004
fp.ReadLine()
# CAR
fp.ReadLine()
# ZnO [PAW_PBE54] (band)
inf["title"] = fp.ReadLine().strip()
# 36 204 24
aa = fp.ReadLine().split()
inf["n"] = pint(aa[0])
inf["nk"] = pint(aa[1])
inf["nLevels"] = pint(aa[2])
Elist = []
Elist_scf = []
Elist_all = []
kx0 = None
ky0 = None
kz0 = None
for ik in range(inf["nk"]):
# next lines can be brank, to be skipped
while 1:
pos = fp.Tell()
line = fp.ReadLine()
if line == '':
break
# print("line1=", line)
if line.strip() != '':
fp.Seek(pos)
break
#
# fp.ReadLine()
# 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.4901961E-02
line = fp.ReadLine()
# print("line2=", line)
if not line: return inf
aa = line.split()
if len(aa) < 4: return inf
kx = pfloat(aa[0])
ky = pfloat(aa[1])
kz = pfloat(aa[2])
wk = pfloat(aa[3])
# 1 -16.924278 1.000000
Eups = []
Edns = []
occups = []
occdns = []
for iL in range(inf["nLevels"]):
line = fp.ReadLine()
# print("line2=", ISPIN, line)
aa = line.split()
aa[1] = pfloat(aa[1])
if len(aa) > 2:
aa[2] = pfloat(aa[2])
else:
if ISPIN == 1:
aa.append(0.0)
if ISPIN == 1 or len(aa) <= 3:
if normalize_E:
Eups.append(aa[1] - EF)
else:
Eups.append(aa[1])
occups.append(aa[2])
else:
if normalize_E:
Eups.append(aa[1] - EF)
Edns.append(aa[2] - EF)
else:
Eups.append(aa[1])
Edns.append(aa[2])
occups.append(pfloat(aa[3]))
occdns.append(pfloat(aa[4]))
dk = None
ktot = None
kinf = [kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns]
Elist_all.append(kinf)
if IsHF:
if not kptinf or kptinf[ik]['kw'] == 0.0:
Elist.append(kinf)
else:
Elist_scf.append(kinf)
else:
Elist.append(kinf)
fp.Close()
# Calculate k distances
ktot = 0.0
for ik in range(len(Elist)):
kx, ky, kz, wk, _dk, _ktot, Eups, occups, Edns, occdns = Elist[ik]
if kx0 is None:
kx0, ky0, kz0 = kx, ky, kz
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx, ky, kz])
ktot += dk
kx0 = kx
ky0 = ky
kz0 = kz
Elist[ik] = [kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns]
ktot = 0.0
for ik in range(len(Elist_scf)):
kx, ky, kz, wk, _dk, _ktot, Eups, occups, Edns, occdns = Elist_scf[ik]
if kx0 is None:
kx0, ky0, kz0 = kx, ky, kz
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx, ky, kz])
ktot += dk
kx0 = kx
ky0 = ky
kz0 = kz
Elist_scf[ik] = [kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns]
ktot = 0.0
for ik in range(len(Elist_all)):
kx, ky, kz, wk, _dk, _ktot, Eups, occups, Edns, occdns = Elist_all[ik]
if kx0 is None:
kx0, ky0, kz0 = kx, ky, kz
dk = cry.GetReciprocalDistanceFromK([kx0, ky0, kz0], [kx, ky, kz])
ktot += dk
kx0 = kx
ky0 = ky
kz0 = kz
Elist_all[ik] = [kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns]
inf['EList'] = Elist
inf['EList_scf'] = Elist_scf
inf['EList_all'] = Elist_all
self.eigenvalinf = inf
return inf
[ドキュメント]
def find_band_edges_from_dos(self, E, DOS, EF0 = 0.0, DOSth = 1.0e18):
"""
DOS（状態密度）データから価電子帯上端 (EV) と伝導帯下端 (EC) を見つけます。
フェルミエネルギー (EF0) を基準として、指定されたDOSスレッショルド (DOSth) を超える
DOS値を持つエネルギー点を探索し、EVとECを決定します。
:param E: numpy.ndarray: エネルギー点の配列。
:param DOS: numpy.ndarray: 対応する状態密度の配列。
:param EF0: float: フェルミエネルギー。デフォルトは0.0。
:param DOSth: float: DOSの閾値。この値を超えるとバンドが存在すると見なされる。デフォルトは1.0e18。
:returns: tuple: (float, float): 価電子帯上端 (EV) と伝導帯下端 (EC)。
"""
EV = 1.0e30
EC = 1.0e30
i0 = -1
for i in range(len(E)):
# print("i=", i, E[i], EF0)
if E[i] >= EF0:
i0 = i
break
for i in range(i0, 0, -1):
# print("i2=", i, DOS[i], DOSth)
if DOS[i] > DOSth:
EV = E[i+1]
break
for i in range(i0, len(E), +1):
# print("i3=", i, DOS[i], DOSth)
if DOS[i] > DOSth:
EC = E[i-1]
break
if EC < EV:
EC = EF0
EV = EF0
return EV, EC
[ドキュメント]
def gbandedges(self, dir = '.', fn_outcar = 'OUTCAR', fn_eig = 'EIGENVAL', ef = None):
"""
EIGENVALファイルとOUTCARファイルからバンドギャップ情報を抽出します。
EIGENVALファイルから各バンドのエネルギーを読み込み、OUTCARファイルから取得した
フェルミエネルギーを基準にして、価電子帯上端 (EVBM) と伝導帯下端 (ECBM) を特定します。
それらの情報とK点座標、バンドインデックス、スピンなどの詳細を辞書で返します。
:param dir: str: VASP計算ディレクトリのパス。デフォルトは'.'。
:param fn_outcar: str: OUTCARファイル名。デフォルトは'OUTCAR'。
:param fn_eig: str: EIGENVALファイル名。デフォルトは'EIGENVAL'。
:param ef: float or None: フェルミエネルギー。Noneの場合、OUTCARから読み込む。
:returns: dict or None: バンドギャップ情報を含む辞書、またはEIGENVALファイルが見つからない場合はNone。
"""
inf = {}
"""
if ef is None:
ief = 0
else:
ief = 1
if not ief:
if os.path.exists(fn_outcar):
for line in open(fn_outcar, 'r'):
if line.find('E-fermi') > -1:
ef = float(line.split()[2])
else:
print(f'Error in tkVASP.gbandedges: [{fn_outcar}] does not exist!')
return None
"""
if not os.path.exists(fn_eig):
print(f'Error in tkVASP.gbandedges: [{fn_eig}] does not exist!')
return None
f_eig = open(fn_eig)
lines = f_eig.readlines()
f_eig.close()
nspin = int(lines[0].split()[3])
nkpoint = int(lines[5].split()[1])
nband = int(lines[5].split()[2])
eup = ef + 99.0
edown = ef - 99.0
k = []
for i in range(nkpoint):
iline = 7 + i * (nband + 2)
# print(f"lines[{iline}]: {lines[iline]}", end = '')
k.append(lines[iline].split()[0:3])
for j in range(nband):
for s in range(nspin):
eval = float(lines[8 + i * (nband + 2) + j].split()[s + 1])
d = eval - ef
if (d > 0):
if (eval < eup):
eup = eval
kup = i
bup = j
sup = s
if (d < 0):
if (eval > edown):
edown = eval
kdown = i
bdown = j
sdown = s
inf["nspin"] = nspin
inf["nkpoint"] = nkpoint
inf["nband"] = nband
inf["nkpoint"] = nkpoint
inf["nlines"] = len(lines)
inf["EF"] = ef
inf["Eg"] = eup - edown
inf["ECBM"] = eup
inf["spin_CBM"] = sup + 1
inf["ik_CBM"] = kup+1
inf["k_CBM"] = [float(k) for k in k[kup]]
inf["iband_CBM"] = bup + 11
inf["spin_CBM"] = sup + 1
inf["EVBM"] = edown
inf["ik_VBM"] = kdown + 1
inf["k_VBM"] = [float(k) for k in k[kdown]]
inf["iband_VBM"] = bdown + 11
inf["spin_VBM"] = sdown + 1
# print("inf=", inf)
return inf
[ドキュメント]
def find_band_edges_from_eigenval(self, EF0 = 0.0, normalize_E = True, eigenvalinf = None, ISPIN = None,
occ_th = 0.5, occ_th0 = 1.0e-4, print_level = 1):
"""
バンドエネルギー情報から価電子帯上端 (EV)、伝導帯下端 (EC)、HOMO、LUMOを見つけます。
EIGENVALから読み込んだバンドエネルギーと占有数を用いて、
電子占有の閾値 (`occ_th`, `occ_th0`) を基準に、各K点でのバンド端エネルギーを探索し、
バンドギャップ関連の情報を辞書として返します。
:param EF0: float: フェルミエネルギー。デフォルトは0.0。
:param normalize_E: bool: エネルギーがフェルミ準位で正規化されているか。デフォルトはTrue。
:param eigenvalinf: dict or None: EIGENVAL情報を含む辞書。Noneの場合、オブジェクトの内部情報を使用。
:param ISPIN: int or None: スピン偏極計算 (`ISPIN=2`) かどうか。Noneの場合、OUTCARまたはINCAR情報から取得。
:param occ_th: float: 占有数の閾値（VBM/CBM検出用）。デフォルトは0.5。
:param occ_th0: float: 占有数の閾値（HOMO/LUMO検出用、occ_thより小さい値）。デフォルトは1.0e-4。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict: バンド端情報（EV, EC, Eg, EF0, HOMO, LUMO）を含む辞書。
"""
if ISPIN is None:
if self.outcarinf:
ISPIN = self.outcarinf["ISPIN"]
elif self.incarinf:
ISPIN = self.incarinf["ISPIN"]
else:
ISPIN = 1
if EF0 is None:
EF0 = 0.0
if eigenvalinf is None:
eigenvalinf = self.eigenvalinf
EPS = 1.0e-3
nk = eigenvalinf["nk"]
nLevels = eigenvalinf["nLevels"]
# print("k points in EIGENVAL:")
# print("nk=", nk)
# print("nLevels=", nLevels)
EV = -1e10
EC = 1e10
EHOMO = -1.0e10
ELUMO = 1.0e10
EList = eigenvalinf.get('EList', None)
if EList is None: return {}
for i in range(len(EList)):
el = EList[i]
kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el
if print_level:
print(f" ({kx:8.4f} {ky:8.4f} {kz:8.4f}) w={wk:8.4g} {dk:8.4f} {ktot:12.4f}")
for il in range(nLevels):
if occups[il] > occ_th0 and EHOMO < Eups[il]:
EHOMO = Eups[il]
if occups[il] > occ_th and EV < Eups[il]:
# if Eups[il] < EF0 + EPS and EV < Eups[il]:
EV = Eups[il]
if occups[il] <= occ_th0 and ELUMO > Eups[il]:
ELUMO = Eups[il]
if occups[il] <= occ_th and EC > Eups[il]:
# if Eups[il] > EF0 - EPS and EC > Eups[il]:
EC = Eups[il]
# print(f" il={il} E={Eups[il]} occ={occups[il]} EF0={EF0}")
# print(" CBM")
if ISPIN == 2 and len(occdns) > 1:
if occdns[il] > occ_th0 and EHOMO < Edns[il]:
EHOMO = Edns[il]
if occups[il] > occ_th and EV < Edns[il]:
# if Edns[il] < EF0 + EPS and EV < Edns[il]:
EV = Edns[il]
if occdns[il] <= occ_th0 and ELUMO > Edns[il]:
ELUMO = Edns[il]
if occups[il] <= occ_th and EC > Edns[il]:
# if Edns[il] > EF0 - EPS and EC > Edns[il]:
EC = Edns[il]
# print("edges:", EHOMO, ELUMO, EV, EC)
if EC < EV:
EC = EF0
EV = EF0
EHOMO = EF0
ELUMO = EF0
# print("edges:", EHOMO, ELUMO, EV, EC)
inf = { "EV": EV, "EC": EC, "Eg": EC - EV, "EF0": EF0, "EHOMO": EHOMO, "ELUMO": ELUMO }
return inf
# lm: 'tot'|'up'|'dn': Total DOS, atom_type + 'tot'|'up'|'dn': Atom-decomposed DOS
# 's' 'p', 'd' + ''|'up'|'dn': l-decomposed DOS 's', 'px', etc + ''|'up'|'dn': lm-decomposed DOS
[ドキュメント]
def cal_lmdecomposed_dos(self, inf, atom_type, lm, save_lmdos = False, outfile = None, is_print = False):
"""
DOSCAR情報から特定の原子種と軌道量子数 (lm) の部分状態密度 (PDOS) を計算します。
`read_doscar` で読み込まれたPDOSデータを用いて、指定された `atom_type` と `lm`
（例: 's', 'p', 'd', 's up'など）に対応するDOSを合計します。
結果をファイルに保存するオプションもあります。
:param inf: dict: `read_doscar` から得られたDOSCAR情報を含む辞書。
:param atom_type: str: PDOSを計算する原子のタイプ（例: 'O', 'Fe'）。
:param lm: str: PDOSを計算する軌道量子数またはスピン成分（例: 's', 'p', 'd', 'tot', 'up', 'dn', 's up'）。
:param save_lmdos: bool: 結果をファイルに保存するかどうか。デフォルトはFalse。
:param outfile: str or None: 結果を保存するファイル名。Noneの場合、デフォルト名が使用される。
:param is_print: bool: 計算過程を出力するかどうか。デフォルトはFalse。
:returns: numpy.ndarray or None: 計算されたlm分解DOSの配列、またはデータが存在しない場合はNone。
"""
npdos = inf["npdos"]
atom_types = inf["atom_types"]
atom_sites = inf["atom_sites"]
lmlabels = inf["lmlabels"]
ks = inf["ks"]
ndata = inf["ndata"]
#pdos_atom_list[iatomsite][lm][idata]
pdos_atom_list = inf["PDOS_atom_list"]
if is_print:
print(f"Calculate lm-deomposed DOS for atom [{atom_type}] lm=[{lm}]")
if npdos == 0:
if is_print:
print(f" No lm-decomposed data (npdos={npdos}. Skip")
return None
lmdos = np.zeros(ndata)
for ia in range(len(atom_sites)):
atom_name = atom_sites[ia]
skip = True
choise_total = ['tot', 'up', 'dn']
if lm in choise_total:
skip = False
choise_atom = [f'{atom_name} tot', f'{atom_name} up', f'{atom_name} dn']
if lm in choise_atom or atom_name == atom_type:
skip = False
if skip:
continue
for iE in range(ndata):
for ilm in range(len(lmlabels)):
lmstr = lmlabels[ilm]
# print(f"at={atom_type}/{atom_name} ilm={ilm} lm=[{lm}]/[{lmstr}]")
skip = True
# total/up/dnを通す
if lm == 'tot' or lm == f'{atom_name} tot':
skip = False
elif (lm == 'up' or lm == f'{atom_name} up') and 'up' in lmstr:
skip = False
elif (lm == 'dn' or lm == f'{atom_name} dn') and 'dn' in lmstr:
skip = False
# lm=s up, px up etcを通す
if lm == lmstr:
skip = False
# lm=s, p, d, px, py, etcを通す => s up, s dn, sを通す
elif lm in lmstr:
skip = False
if skip:
continue
# print(" pass")
if lm == 'tot up' and 'up' in lmstr:
lmdos[iE] += ks * pdos_atom_list[ia][ilm][iE]
elif lm == 'tot dn' and 'dn' in lmstr:
lmdos[iE] += ks * pdos_atom_list[ia][ilm][iE]
else:
lmdos[iE] += ks * pdos_atom_list[ia][ilm][iE]
if save_lmdos:
# print("lmdos=", lmdos)
data_list = inf["data_list"][0]
xE = data_list[0]
pdos_atom_list[ia][ilm][iE]
Elabels = ["E (eV)", f"{atom_type} {lm}"]
ylist = [xE, lmdos]
for ia in range(len(atom_sites)):
atom_name = atom_sites[ia]
if lm != 'tot' and atom_name != atom_type:
continue
for ilm in range(len(lmlabels)):
lmstr = lmlabels[ilm]
if lm != 'tot' and lm != f'{atom_type} tot' and lm not in lmstr:
continue
Elabels.append(f"{atom_type}#{ia} {lmstr}#{ilm}")
ylist.append(pdos_atom_list[ia][ilm])
if outfile is None:
if lm == 'tot':
outfile = f"pdos_total.xlsx"
else:
outfile = f"pdos_{atom_type}_{lm}.xlsx"
tkVariousData().to_excel(outfile, Elabels, ylist)
return lmdos
[ドキュメント]
def read_doscar(self, DOSCAR_path, cry = None, EF = None, normalize_E = True, unit = '',
IsSpinPolarized = None, IsNonCollinear = None, SpinOrbit = None, print_level = 1):
"""
DOSCARファイルから状態密度 (DOS) 情報を読み込みます。
DOSCARファイルの内容を解析し、全状態密度 (TDOS) と部分状態密度 (PDOS) を抽出します。
フェルミエネルギーによるエネルギーの正規化、単位変換、スピン偏極や非共線計算の考慮、
PDOSの軌道分解情報などを処理します。
:param DOSCAR_path: str: DOSCARファイルへのパス。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。Noneの場合、内部でPOSCARから読み込む。
:param EF: float or None: フェルミエネルギー。Noneの場合、OUTCAR情報から取得。
:param normalize_E: bool: エネルギーをフェルミ準位で正規化するかどうか。デフォルトはTrue。
:param unit: str: DOSの単位。'/cm3'を指定すると体積で割って正規化。デフォルトは''（無単位）。
:param IsSpinPolarized: bool or None: スピン偏極計算かどうか。Noneの場合、OUTCAR情報から取得。
:param IsNonCollinear: bool or None: 非共線計算かどうか。Noneの場合、OUTCAR情報から取得。
:param SpinOrbit: bool or None: スピン軌道結合計算かどうか。Noneの場合、OUTCAR情報から取得。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict: DOS情報を含む辞書。
"""
# if self.outcarinf is None: return None
if EF is None:
if self.outcarinf:
EF = self.outcarinf.get("EF", 0.0)
else:
EF = 0.0
if cry is None:
cry = self.crystal
ISPIN = 1
if self.outcarinf:
if IsSpinPolarized is None:
ISPIN = self.outcarinf.get("ISPIN", 1)
if ISPIN == 1:
IsSpinPolarized = 0
else:
IsSpinPolarized = 1
if SpinOrbit is None:
LSORBIT = self.outcarinf.get("LSORBIT", 0)
if 'T' in LSORBIT:
SpinOrbit = 1
else:
SpinOrbit = 0
if IsNonCollinear is None:
LNONCOLLINEAR = self.outcarinf.get("LNONCOLLINEAR", 0)
if 'T' in LNONCOLLINEAR:
IsNonCollinear = 1
else:
IsNonCollinear = 0
else:
IsSpinPolarized = 0
LSORBIT = 'F'
SpinOrbit = 0
LNONCOLLINEAR = 'F'
IsNonCollinear = 0
poscar_path = self.get_POSCAR(DOSCAR_path)
if print_level:
print(f"Read crystal structure from [{poscar_path}]")
poscar_inf = self.read_poscar_inf(poscar_path, print_level = 0)
cry = poscar_inf["crystal"]
atom_types = poscar_inf["AtomTypes"]
sites = cry.ExpandedAtomSiteList()
atom_sites = []
for site in sites:
atom_sites.append(site.AtomNameOnly())
natom_types = len(atom_types)
natom_sites = len(atom_sites)
if print_level:
print("natom_types=", natom_types)
print("atom_types =", atom_types)
print("natom_types=", natom_sites)
print("atom_sites =", atom_sites)
print("natom_sites=", natom_sites)
inf = {}
if print_level:
print()
print(f"Read DOSCAR from [{DOSCAR_path}]")
fp = tkFile(DOSCAR_path, 'r')
if not fp or fp.fp is None:
# terminate("Error in read_average_core_potentials: Can not read [{}]".format(DOSCAR_path))
return { "nE": 0, "E": [], "TotalDOSup": [], "TotalDOSdn": [], "Neup": [], "Nedn": [] }
aa = fp.readline().split()
nions_with_empty = pint(aa[0])
nions = pint(aa[1])
include_PDOS = pint(aa[2])
inf["NCDIJ"] = pint(aa[3])
if print_level:
print("nions=", nions)
print("nions_with_empty=", nions_with_empty)
if nions != natom_sites:
print()
print(f" ERROR: nions={nions} taken from {DOSCAR_path} is different from "
+ f"natom_sites={natom_sites} taken from {poscar_path}")
input("Press ENTER to terminate>>\n")
exit()
aa = fp.readline().split()
Vcell = pfloat(aa[0])
a, b, c = pfloat(aa[1]), pfloat(aa[2]), pfloat(aa[3])
inf["POTIM"] = pfloat(aa[4])
line = fp.readline()
TEBEG = pfloat(line)
fp.readline()
line = fp.readline()
sample_name = line.strip()
# print("s=", inf["SampleName"])
# t/PDOSの読み込み
if print_level >= 2:
print("Read total DOS")
line = fp.ReadLine()
aa = line.split()
if len(aa) < 5:
#0123456789012345678901230123456789012345678901234567890123456789012345678901234567890123456789
# 10.00000000 -30.0000000040001 -0.73600000 1.00000000
# print(f"lines=[{line[0:16]}][{line[16:32]}][{line[32:37]}][{line[37:53]}][{line[53:]}]")
Emax = pfloat(line[0:16])
Emin = pfloat(line[16:32])
nE = pint(line[32:37])
Efermi = pfloat(line[37:53])
# print("E=", Emax, Emin, nE, Efermi)
else:
Emax = pfloat(aa[0])
Emin = pfloat(aa[1])
nE = pint(aa[2])
Efermi = pfloat(aa[3])
# 指数が3桁の場合、eが削除されるので追加
def vfloat(s):
m = re.sub(r"(\d)-(\d)", "$1e-$2$", s)
if m:
s = m
return pfloat(s)
ylist_total = None
for i in range(nE):
line = fp.readline()
if not line:
break
aa = line.split()
if len(aa) == 0:
break
nc = len(aa)
if ylist_total is None:
ylist_total = np.zeros([nc, nE])
for ic in range(nc):
ylist_total[ic][i] = vfloat(aa[ic])
# print("e=", ylist[0][i], ylist[1][i], ylist[2][i])
if ylist_total is None: return inf
# lm-decomposed DOSの読み込み
if print_level >= 2:
print("Read lm-decomposed DOS")
nPDOS = 0
ylist_lmdos = []
while 1:
line = fp.ReadLine()
if not line:
break
if print_level >= 2:
print(f" Read lm-decomposed DOS for {nPDOS}-th site ({atom_sites[nPDOS]})")
ylist = None
for idata in range(nE):
line = fp.readline()
if not line:
break
aa = line.split()
if len(aa) == 0:
break
ncolumn = len(aa)
if ylist is None:
ylist = np.zeros([ncolumn - 1, nE])
for ic in range(1, ncolumn):
ylist[ic - 1][idata] = vfloat(aa[ic])
ylist_lmdos.append(ylist)
nPDOS += 1
# SpinOrbit == 1 の場合、さらに mx, my, mzのデータが続く
fp.Close()
if len(ylist_lmdos) == 0 or ylist_lmdos[0] is None:
ylist_lmdos = [[], []]
npdos = 0
else:
npdos = len(ylist_lmdos[0])
if normalize_E:
Elist = [e - EF for e in ylist_total[0]]
else:
Elist = ylist_total[0]
Vcell = cry.Volume()
k = 1.0 / (Vcell * 1.0e-24)
if unit == '/cm3':
for i in range(1, len(ylist_total)):
ylist_total[i] *= k
if npdos > 0:
for isite in range(nions):
for ilm in range(npdos):
# print("i=", isite, ilm)
ylist_lmdos[isite][ilm] *= k
if nc <= 3:
tDOSuplist = ylist_total[1]
tDOSdnlist = None
Neuplist = ylist_total[2]
Nednlist = None
else:
tDOSuplist = ylist_total[1]
tDOSdnlist = ylist_total[2]
Neuplist = ylist_total[3]
Nednlist = ylist_total[4]
PDOS_atom_list = ylist_lmdos
# total / sum DOSの計算
if tDOSdnlist is None:
tDOSlist = tDOSuplist
Nelist = Neuplist
else:
tDOSlist = np.zeros(nE)
Nelist = np.zeros(nE)
for i in range(nE):
tDOSlist[i] = tDOSuplist[i] + tDOSdnlist[i]
Nelist[i] = Neuplist[i] + Nednlist[i]
if tDOSdnlist is None:
if npdos == 3:
lmlabels = ["s", "p", "d"]
else:
lmlabels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxy", "dx2-y2"]
else:
if npdos == 6:
lmlabels = ["s up", "s dn", "p up", "p dn", "d up", "d dn"]
else:
lmlabels = ["s up", "s dn",
"py up", "py dn", "pz up", "pz dn", "px up", "px dn",
"dxy up", "dxy dn", "dyz up", "dyz dn", "dz2 up", "dz2 dn", "dxz up", "dxz dn", "dx2-y2 up", "dx2-y2 dn"]
# 原子種ごとにPDOS (lmlabels全て) の和を取る
# IsSpinPolarized == 0の場合、2倍する
if SpinOrbit == 1 and IsSpinPolarized == 0:
ks = 2.0
else:
ks = 1.0
PDOS_list = {}
sumDOS_list = np.zeros(nE)
if npdos > 0:
for iion in range(nions):
atom_name = atom_sites[iion]
PDOS_list[atom_name] = np.zeros(nE)
for iion in range(nions):
atom_name = atom_sites[iion]
for ilm in range(len(lmlabels)):
if len(PDOS_atom_list[iion]) <= ilm: continue
for iE in range(nE):
v = ks * PDOS_atom_list[iion][ilm][iE]
PDOS_list[atom_name][iE] += v
sumDOS_list[iE] += v
interstitialDOS_list = tDOSlist - sumDOS_list
else:
interstitialDOS_list = np.zeros(nE)
inf["cry"] = cry
inf["SampleName"] = sample_name
inf["IsSpinPolarized"] = IsSpinPolarized
inf["ISPIN"] = ISPIN
inf["SpinOrbit"] = SpinOrbit
inf["IsNonCollinear"] = IsNonCollinear
inf["ks"] = ks
inf["nions_with_empty"] = nions_with_empty
inf["nions"] = nions
inf["atom_types"] = atom_types
inf["atom_sites"] = atom_sites
inf["natom_types"] = natom_types
inf["natom_sites"] = natom_sites
inf["Vcell"] = Vcell
inf["a"] = a
inf["b"] = b
inf["c"] = c
inf["EF"] = EF
inf["Efermi"] = Efermi
inf["TEBEG"] = TEBEG
inf["include_PODS"] = include_PDOS
inf["Emax"] = Emax
inf["Emin"] = Emin
inf["nE"] = nE
inf["E"] = Elist
inf["TotalDOS"] = tDOSlist
inf["TotalDOSup"] = tDOSuplist
inf["TotalDOSdn"] = tDOSdnlist
inf["Ne"] = Nelist
inf["Neup"] = Neuplist
inf["Nedn"] = Nednlist
inf["npdos"] = npdos
inf["lmlabels"] = lmlabels
inf["PDOS_list"] = PDOS_list
# inf["PDOS_list_labels"] = atom_types
inf["PDOS_atom_list"] = PDOS_atom_list
inf["sumDOS_list"] = sumDOS_list
include_sum = False
if tDOSdnlist is None:
if include_sum:
data_list = [Elist, tDOSlist, sumDOS_list, interstitialDOS_list]
labels = ["E (eV)", "Total", "sum", "interstitial"]
else:
data_list = [Elist, tDOSlist, interstitialDOS_list]
labels = ["E (eV)", "Total", "interstitial"]
if npdos > 0:
for atom_name in atom_types:
data_list.extend([PDOS_list[atom_name]])
labels.extend(atom_types)
else:
if include_sum:
data_list = [Elist, tDOSlist, tDOSuplist, tDOSdnlist, sumDOS_list, interstitialDOS_list]
labels = ["E (eV)", "Total", "Total up", "Total dn", "sum", "interstitial"]
else:
data_list = [Elist, tDOSlist, tDOSuplist, tDOSdnlist, interstitialDOS_list]
labels = ["E (eV)", "Total", "Total up", "Total dn", "interstitial"]
if npdos > 0:
for atom_name in atom_types:
data_list.extend([PDOS_list[atom_name]])
labels.extend(atom_types)
inf["data_list"] = data_list
inf["labels"] = labels
inf["lm_labels"] = lmlabels
# for id in range(len(data_list)):
# print("id=", id, labels[id], type(data_list[id]))
# for iE in range(len(data_list[id])):
# print(" ", iE, data_list[id][iE])
return inf
[ドキュメント]
def read_DOS(self, file, cry = None, IsSpinPolarized = None, IsNonCollinear = None, EF = None, normalize_E = True, unit = '', print_level = 1):
"""
DOSファイルを読み込み、状態密度情報を取得します。
`.dat` ファイル形式またはDOSCARファイル形式のDOSデータを読み込み、
エネルギーと状態密度の情報を辞書形式で返します。
DOSCARの場合、`read_doscar` メソッドを呼び出します。
:param file: str: DOSファイルへのパス。
:param cry: tkCrystal or None: 結晶構造情報を持つ `tkCrystal` オブジェクト。
:param IsSpinPolarized: bool or None: スピン偏極計算かどうか。
:param IsNonCollinear: bool or None: 非共線計算かどうか。
:param EF: float or None: フェルミエネルギー。
:param normalize_E: bool: エネルギーをフェルミ準位で正規化するかどうか。デフォルトはTrue。
:param unit: str: DOSの単位。デフォルトは''（無単位）。
:param print_level: int: 出力レベル。デフォルトは1。
:returns: dict: DOS情報を含む辞書。
"""
inf = {}
if '.dat' in file:
E_raw, dos_raw = np.genfromtxt(file, skip_header=1, usecols=(0,1), unpack=True)
inf["E"] = E_raw
inf["nE"] = len(E_raw)
inf["TotalDOS"] = dos_raw
else:
inf = self.read_doscar(file, cry = cry, EF = EF, normalize_E = normalize_E, unit = unit,
IsSpinPolarized = IsSpinPolarized, IsNonCollinear = IsNonCollinear, print_level = print_level)
return inf
[ドキュメント]
def read_vasprun(self, VASPRUN_path, print_level = 1):
"""
指定されたvasprun.xmlファイルを読み込み、その内容を辞書型で返します。
`xmltodict` ライブラリを使用して `vasprun.xml` を解析し、XML構造をPythonの辞書に変換します。
数値型や真偽値型のデータは適切な型に自動変換されます。
:param VASPRUN_path: str: vasprun.xmlファイルのパス。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは1。
:returns: dict or None: XMLファイルの内容を表す辞書。
ファイルが見つからない場合やパースエラーが発生した場合はNoneを返します。
"""
def _datatype_postprocessor(path, key, value):
# dict で @datatype と #text を含む場合に型変換
if isinstance(value, dict) and '@datatype' in value and '#text' in value:
dt = value.get('@datatype', '').lower()
txt = value.get('#text')
try:
if dt in ('float', 'double', 'real'):
new_val = float(txt)
elif dt in ('int', 'integer', 'i'):
new_val = int(txt)
elif dt in ('bool', 'boolean'):
new_val = txt.lower() in ('true', '1')
else:
new_val = txt
except Exception:
new_val = txt
return key, new_val
return key, value
# parse_dos や parse_eigen を True にすると DOS/EIGEN 情報も読み込みます
path = Path(VASPRUN_path)
raw_xml = path.read_text(encoding='utf-8')
# xmltodict で解析、postprocessor で型変換
parsed = xmltodict.parse(
raw_xml,
attr_prefix='@',
cdata_key='#text',
postprocessor=_datatype_postprocessor
)
return parsed
"""
with open(VASPRUN_path, 'rb') as f:
xml_bytes = f.read()
d = xmltodict.parse(xml_bytes)
return d
"""
[ドキュメント]
def read_files(self, file, file_types, EF = None, normalize_E = True, unit = '/cm3', data_for_bandedges = None,
exit_by_error = False, print_level = 0, terminate = terminate):
"""
指定されたVASP計算ディレクトリ内の複数のVASP出力ファイルを一括で読み込みます。
INCAR, POSCAR, POTCAR, KPOINTS, CONTCAR, OUTCAR, EIGENVAL, DOSCAR, vasprun.xml などのファイルを
指定された `file_types` に応じて読み込みます。各ファイルの解析結果は、返される辞書に格納されます。
また、計算の役割（scf, band, relaxなど）も推測します。
:param file: str: VASP計算ディレクトリ内にあるファイルまたはディレクトリのパス。このパスから計算ディレクトリを特定します。
:param file_types: list of str: 読み込むファイルの種類を指定する文字列のリスト（例: ['INCAR', 'OUTCAR', 'DOSCAR']）。
:param EF: float or None: フェルミエネルギー。Noneの場合、OUTCAR情報から取得。
:param normalize_E: bool: エネルギーをフェルミ準位で正規化するかどうか。デフォルトはTrue。
:param unit: str: DOSの単位。'/cm3'を指定すると体積で割って正規化。デフォルトは'/cm3'。
:param data_for_bandedges: Any: バンド端計算のための追加データ。現在は使用されていません。
:param exit_by_error: bool: ファイル読み込みエラー時にプログラムを終了するかどうか。デフォルトはFalse。
:param print_level: int: 出力レベル。0で非表示、1以上で表示。デフォルトは0。
:param terminate: callable: エラー時に呼び出される終了関数。デフォルトはtklib.tkutils.terminate。
:returns: dict: 読み込まれたVASPファイルの情報を含む辞書。
"""
base_path = self.getdir(file)
INCAR_path = self.get_INCAR(base_path)
POSCAR_path = self.get_POSCAR(base_path)
POTCAR_path = self.get_POTCAR(base_path)
KPOINTS_path = self.get_KPOINTS(base_path)
CONTCAR_path = self.get_CONTCAR(base_path)
OUTCAR_path = self.get_OUTCAR(base_path)
EIGENVAL_path = self.get_VASPPath(base_path, 'EIGENVAL')
EIGENVAL_OPT_path = self.get_VASPPath(base_path, 'EIGENVAL_OPT')
DOSCAR_path = self.get_VASPPath(base_path, 'DOSCAR')
VASPRUN_path = self.get_VASPRUN(base_path)
if print_level >= 1:
print("")
print("CAR dir : ", base_path)
print(" INCAR : ", INCAR_path)
print(" POSCAR : ", POSCAR_path)
print(" POTCAR : ", POTCAR_path)
print(" KPOINTS : ", KPOINTS_path)
print(" CONTCAR : ", CONTCAR_path)
print(" OUTCAR : ", OUTCAR_path)
print(" EIGENVAL: ", EIGENVAL_path)
print(" EIGENVAL_OPT: ", EIGENVAL_OPT_path)
print(" DOSCAR : ", DOSCAR_path)
print(" vasprun.xml : ", VASPRUN_path)
abs_base_path = os.path.abspath(os.path.dirname(base_path))
last_dir = os.path.basename(abs_base_path.rstrip('/'))
parent_dir = os.path.dirname(abs_base_path)
inf = { "base_path": base_path, "last_dir": last_dir, "parent_dir": parent_dir, "role": ""}
if True:
if print_level >= 1:
print("*** Read INCAR from [{}]".format(INCAR_path))
_inf = self.read_incar_inf(INCAR_path, print_level = print_level)
if exit_by_error and _inf is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(INCAR_path))
if 'INCAR' in file_types:
inf["INCAR"] = _inf
if 'POTCAR' in file_types:
if print_level >= 1:
print("*** Read POTCAR from [{}]".format(POTCAR_path))
inf["POTCAR"] = self.read_potcar_inf(POTCAR_path, print_level = print_level)
if exit_by_error and inf["POTCAR"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(POTCAR_path))
if 'POSCAR' in file_types:
if print_level >= 1:
print("*** Read crystal structure from [{}]".format(POSCAR_path))
inf["POSCAR"] = self.read_poscar_inf(POSCAR_path, print_level = print_level)
if exit_by_error and inf["POSCAR"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(POSCAR_path))
if True:
if print_level >= 1:
print(f"*** Read k points from [{KPOINTS_path}]")
inf_kpoints = self.read_kpoints(KPOINTS_path)
if exit_by_error and inf_kpoints is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(KPOINTS_path))
if 'KPOINTS' in file_types:
inf["KPOINTS"] = inf_kpoints
if 'CONTCAR' in file_types:
if print_level >= 1:
print("*** Read crystal structure from [{}]".format(CONTCAR_path))
inf["CONTCAR"] = self.read_poscar_inf(CONTCAR_path, print_level = print_level)
if exit_by_error and inf["POSCAR"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(CONTCAR_path))
if True:
if print_level >= 1:
print("*** Read OUTCAR from [{}]".format(OUTCAR_path))
inf_outcar = self.read_outcar_inf(OUTCAR_path, print_level = print_level)
if exit_by_error and inf_outcar is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(OUTCAR_path))
if 'OUTCAR' in file_types:
inf["OUTCAR"] = inf_outcar
def speculate_role(inf_outcar, inf_kpoints):
ICHARG = inf_outcar.get("ICHARG", 100000)
IBRION = inf_outcar.get("IBRION", 100000)
ISIF = inf_outcar.get("ISIF", 100000)
LKPOINTS_OPT = 'T' in inf_outcar.get("LKPOINTS_OPT", "").upper()
LPARD = 'T' in inf_outcar.get("LPARD", "").upper()
line_mode = 'L' in inf_kpoints.get("mode", "").upper()
if ICHARG < 10 : scf = True
elif ICHARG < 100: scf = False
if IBRION == -1:
if scf:
return "scf"
elif LKPOINTS_OPT or line_mode:
return "band"
elif LPARD:
return "electron density"
else:
return"dos"
elif IBRION == 0:
if ISIF >= 3:
return"vc md"
else:
return "md"
elif 1 <= IBRION <= 3:
if ISIF >= 3:
return "vc relax"
else:
return "relax"
elif 5 <= IBRION <= 8: return "phonon"
elif IBRION in [40, 44]: return "transition state"
elif IBRION in [11, 12]: return "user defined"
return "unknown"
inf["role"] = speculate_role(inf_outcar, inf_kpoints)
if 'VASPRUN' in file_types:
if print_level >= 1:
print("*** Read vasprun.xml from [{}]".format(VASPRUN_path))
inf["VASPRUN"] = self.read_vasprun(VASPRUN_path, print_level = print_level)
if exit_by_error and inf["VASPRUN"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(VASPRUN_path))
if 'DOSCAR' in file_types:
if print_level >= 1:
print("*** Read DOSCAR from [{}]".format(DOSCAR_path))
inf["DOSCAR"] = self.read_doscar(DOSCAR_path, EF = EF, normalize_E = normalize_E, unit = unit, print_level = print_level)
if exit_by_error and inf["DOSCAR"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(DOSCAR_path))
if 'EIGENVAL' in file_types:
if print_level >= 1:
print("*** Read EIGENVAL from [{}]".format(EIGENVAL_path))
inf["EIGENVAL"] = self.read_eigenval(EIGENVAL_path, normalize_E = normalize_E, print_level = print_level)
if exit_by_error and inf["EIGENVAL"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(EIGENVAL_path))
if 'EIGENVAL_OPT' in file_types:
if print_level >= 1:
print("*** Read EIGENVAL_OPT from [{}]".format(EIGENVAL_OPT_path))
inf["EIGENVAL_OPT"] = self.read_eigenval_opt(EIGENVAL_OPT_path, normalize_E = normalize_E, print_level = print_level)
if exit_by_error and inf["EIGENVAL_OPT"] is None:
terminate("Error in tkvasp.read_files: Can not read [{}]".format(EIGENVAL_OPT_path))
return inf