"""
VASPの計算結果から帯電欠陥および双極子相互作用の補正を推定するスクリプト。
このスクリプトは、VASPの計算結果（OUTCARやINCARなど）を解析し、
帯電欠陥を持つ系における長距離クーロン相互作用や双極子相互作用によって生じる
有限サイズ補正エネルギーを推定します。
具体的には、IDIPOL=4の計算結果、電子誘電率、イオン性誘電率、
そして欠陥の電荷情報を用いて補正値を算出します。
使用方法の詳細は以下のドキュメントを参照してください。
:doc:`vasp_correction_charge_usage`
"""
import os
import sys
import shutil
import glob
import csv
import re
import numpy as np
from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, pi
from scipy.interpolate import interp1d
from pprint import pprint
from matplotlib import pyplot as plt
from tklib.tkfile import tkFile
from tklib.tkinifile import tkIniFile
import tklib.tkre as tkre
from tklib.tkapplication import tkApplication
from tklib.tkutils import save_csv
from tklib.tkutils import IsDir, IsFile, SplitFilePath
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tksci.tksci import Reduce01, Round
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
from tklib.tkcrystal.tkcif import tkCIF, tkCIFData
from tklib.tkcrystal.tkcrystal import tkCrystal
from tklib.tkcrystal.tkatomtype import tkAtomType
from tklib.tkcrystal.tkvasp import tkVASP
from tklib.tkcrystal.tkbandstructure import plot_band
from tklib.tksci import tkequation
import tklib.tkcsv
"""
Estimate the correction of charged defects and dipole interaction by VASP
"""
#================================
# global parameters
#================================
debug = 0
validmodes = ['Z']
mode = 'Z'
CAR_IDIPOL_dir = '.'
CAR_DIEL_el_dir = '.'
CAR_DIEL_ion_dir = '.'
CAR_charge_dir = '.'
summaryfile = None
# Energy to search edge energies
EF0 = 0.1
# Max number of electrons occupies one state
Nemax = 1.0
#===================================
# figure configuration
#===================================
band_marker_size = 4
band_marker_edge_width = 0.5
Emin = -10.0 # eV
Emax = 5.0 # eV
figsize = (6, 8)
#figsize = (6, 4)
#colors = ['#000000', '#ff0000', '#00aa00', '#0000ff', '#aaaa00', '#ff00ff', '#00ffff', '#aa0000', '#00aa00', '#0000aa']
#colors = ['k', 'r', 'g', 'b', 'y', 'm', 'c']
fontsize = 16
labelfontsize = 12
legend_fontsize = 8
app = tkApplication()
#=============================
# Treat argments
#=============================
[ドキュメント]
def usage(app = None):
"""
スクリプトの正しい使用方法を標準出力に表示する。
詳細説明:
本スクリプトのコマンドライン引数の渡し方を示す。
引数には、IDIPOL=4計算ディレクトリ、電子誘電率計算ディレクトリまたは値、
イオン誘電率計算ディレクトリまたは値、欠陥電荷計算ディレクトリを渡す。
:param app: tkApplicationのインスタンス（オプション）。
:returns: None
"""
print("")
print("Usage:")
print(" (a) python {} CAR_IDIPOL_dir CAR_DIEL_el_dir CAR_DIEL_ion_dir CAR_charge_dir".format(sys.argv[0]))
print(" (b) python {} CAR_IDIPOL_dir CAR_DIEL_el_dir epsilon CAR_charge_dir".format(sys.argv[0]))
print(" ex: python {} {} {} {} {}".format(sys.argv[0], CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir))
[ドキュメント]
def updatevars():
"""
コマンドライン引数に基づいてグローバル変数を更新する。
詳細説明:
`sys.argv`からコマンドライン引数を取得し、
`CAR_IDIPOL_dir`, `CAR_DIEL_el_dir`, `CAR_DIEL_ion_dir`, `CAR_charge_dir`
のグローバル変数を更新します。
現在サポートされている`mode`は'Z'のみであり、それ以外の値が指定された場合は
エラーメッセージを表示してプログラムを終了します。
:returns: None
"""
global mode
global CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir
# mode = getarg (1, mode)
CAR_IDIPOL_dir = getarg(1, CAR_IDIPOL_dir)
CAR_DIEL_el_dir = getarg(2, CAR_DIEL_el_dir)
CAR_DIEL_ion_dir = getarg(3, CAR_DIEL_ion_dir)
CAR_charge_dir = getarg(4, CAR_charge_dir)
if mode == 'Z':
pass
else:
app.terminate(f"Error: Invalide mode [{mode}]", usage = usage, pause = True)
#def save_csv(path, headerlist, datalist, is_print = 0):
[ドキュメント]
def read_dipole_inf(outcarpath):
"""
VASPのOUTCARファイルから双極子補正に関する情報を読み込む。
詳細説明:
指定されたOUTCARファイルを解析し、"DIPCOR: dipole corrections for dipol"
セクションから以下の情報を抽出します。
- 双極子モーメント (dipolemoment)
- 四極子モーメントのトレース (Tr[quadrupol])
- 電荷によるエネルギー補正 (Ecorr(charge))
- 双極子+四極子によるエネルギー補正 (Ecorr(di+quadrupol))
- 追加電場によるイオン相互作用のエネルギー補正 (Ecorr(with field))
情報が見つからない場合やファイルが読み込めない場合はNoneを返します。
:param outcarpath: str OUTCARファイルのパス。
:returns: dict 読み込んだ双極子補正情報を含む辞書。ファイルが読み込めない場合はNone。
"""
fp = tkFile(outcarpath, 'r')
if not fp or fp.fp is None:
return None
inf = {}
line = fp.SkipTo('DIPCOR: dipole corrections for dipol')
if line:
# direction 1 min pos 1, direction 2 min pos 1, direction 3 min pos 7,
fp.ReadLine()
# dipolmoment 0.000000 -0.000000 -6.034717 electrons x Angstroem
# Tr[quadrupol] -5.863427
line = fp.ReadLine()
# print("l=", line)
aa = line.split()
inf["dipolemoment"] = [pfloat(aa[1]), pfloat(aa[2]), pfloat(aa[3])]
aa = fp.ReadLine().split()
inf["Tr[quadrupol]"] = pfloat(aa[1])
#
fp.ReadLine()
# energy correction for charged system 0.000000 eV
line = fp.ReadLine()
s = tkre.Search(r'([+\-]?[+\-\d\.eEdD]+)\s*eV', line)
inf["Ecorr(charge)"] = pfloat(s[1])
# dipol+quadrupol energy correction -14.435503 eV
line = fp.ReadLine()
s = tkre.Search(r'([+\-]?[+\-\d\.eEdD]+)\s*eV', line)
# print("s=", s)
inf["Ecorr(di+quadrupol)"] = pfloat(s[1])
# added-field ion interaction 0.000000 eV (added to PSCEN)
line = fp.ReadLine()
s = tkre.Search(r'([+\-]?[+\-\d\.eEdD]+)\s*eV', line)
inf["Ecorr(with field)"] = pfloat(s[1])
fp.Close()
return inf
[ドキュメント]
def read_dielectric_inf_ion(outcarpath):
"""
VASPのOUTCARファイルからイオン性誘電率に関する情報を読み込む。
詳細説明:
指定されたOUTCARファイルを解析し、以下のセクションから誘電率テンソルを抽出します。
- "MACROSCOPIC STATIC DIELECTRIC TENSOR (incl"：マクロな静的誘電率（局所電場を含む）
- "MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC"：イオン性誘電率
これらのテンソルを辞書に格納して返します。情報が見つからない場合や
ファイルが読み込めない場合はNoneを返します。
:param outcarpath: str OUTCARファイルのパス。
:returns: dict 読み込んだイオン性誘電率情報を含む辞書。ファイルが読み込めない場合はNone。
"""
fp = tkFile(outcarpath, 'r')
if not fp or fp.fp is None:
return None
inf = {}
while 1:
eps1 = make_matrix2(3, 3)
eps2 = make_matrix2(3, 3)
line = fp.SkipTo(r"MACROSCOPIC STATIC DIELECTRIC TENSOR \(incl")
if not line:
break
# ------------------------------------------------------
# 7.930930 0.000000 0.000000
# 0.000000 7.930923 0.000000
# -0.000000 0.000000 7.351176
# ------------------------------------------------------
line = fp.ReadLine()
aa = make_matrix1(3)
aa[0] = fp.ReadLine().split()
aa[1] = fp.ReadLine().split()
aa[2] = fp.ReadLine().split()
for i in range(3):
eps1[i][0] = pfloat(aa[i][0])
eps1[i][1] = pfloat(aa[i][1])
eps1[i][2] = pfloat(aa[i][2])
inf["epsilon(static/local field)"] = eps1
line = fp.SkipTo("MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC")
if not line:
break
line = fp.ReadLine()
aa[0] = fp.ReadLine().split()
aa[1] = fp.ReadLine().split()
aa[2] = fp.ReadLine().split()
for i in range(3):
eps2[i][0] = pfloat(aa[i][0])
eps2[i][1] = pfloat(aa[i][1])
eps2[i][2] = pfloat(aa[i][2])
inf["epsilon(static/ion)"] = eps2
return inf
[ドキュメント]
def read_dielectric_inf_el(outcarpath):
"""
VASPのOUTCARファイルから電子誘電率に関する情報を読み込む。
詳細説明:
指定されたOUTCARファイルを解析し、"REAL DIELECTRIC FUNCTION"セクションから
実部の誘電率テンソルを抽出します。
通常は、エネルギーが0.000000 eVに対応する最初の誘電率データを使用します。
情報が見つからない場合やファイルが読み込めない場合はNoneを返します。
:param outcarpath: str OUTCARファイルのパス。
:returns: dict 読み込んだ電子誘電率情報を含む辞書。ファイルが読み込めない場合はNone。
"""
fp = tkFile(outcarpath, 'r')
if not fp or fp.fp is None:
return None
inf = {}
while 1:
line = fp.SkipTo(r"REAL DIELECTRIC FUNCTION")
if not line:
break
# E(ev) X Y Z XY YZ ZX
# --------------------------------------------------------------------------------------------------------------
# 0.000000 2.861856 2.864431 2.893541 0.000807 0.000000 0.000000
# 0.058843 2.861930 2.864506 2.893617 0.000807 0.000000 0.000000
# 0.117686 2.862151 2.864728 2.893843 0.000808 0.000000 0.000000
fp.ReadLine()
fp.ReadLine()
line = fp.ReadLine()
aa = line.split()
eps_r = make_matrix2(3, 3)
eps_r[0][0] = pfloat(aa[1]) # X
eps_r[1][1] = pfloat(aa[2]) # Y
eps_r[2][2] = pfloat(aa[3]) # Z
eps_r[0][1] = pfloat(aa[4]) # XY
eps_r[1][0] = pfloat(aa[4]) # XY
eps_r[0][2] = pfloat(aa[6]) # ZX
eps_r[2][0] = pfloat(aa[6]) # ZX
eps_r[1][2] = pfloat(aa[5]) # YZ
eps_r[2][1] = pfloat(aa[5]) # YZ
inf["eps_opt_r"] = eps_r
break
return inf
[ドキュメント]
def correct_charged_defects(mode, CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir):
"""
帯電欠陥および双極子相互作用のエネルギー補正を計算し、結果を要約ファイルに保存する。
詳細説明:
以下の4つのVASP計算ディレクトリから必要な情報を収集し、帯電欠陥のエネルギー補正を計算します。
1. `CAR_IDIPOL_dir`: IDIPOL=4で計算された欠陥モデルのディレクトリ。
ここから双極子補正情報（Ecorr(charge)など）を読み取ります。
2. `CAR_DIEL_el_dir`: 電子誘電率計算モデルのディレクトリ、または電子誘電率の平均値。
電子誘電率を読み取ります。
3. `CAR_DIEL_ion_dir`: イオン性誘電率計算モデルのディレクトリ、またはイオン性誘電率の平均値。
イオン性誘電率を読み取ります。
4. `CAR_charge_dir`: 欠陥の電荷を決定するための欠陥モデルのディレクトリ。
POSCARとPOTCARから総価電子数を取得し、OUTCARからNELECTを読み取って、欠陥の電荷 `q` を決定します。
また、セル体積から欠陥の代表長 `L` を算出します。
集めた情報 (`q`, `L`, `eps_tot`, `dipoleinf["Ecorr(charge)"]`) を使用して、
電荷補正エネルギー `dEcorr` を算出し、結果を`summaryfile`に保存します。
:param mode: str 計算モード（現在'Z'のみ有効）。
:param CAR_IDIPOL_dir: str IDIPOL=4で計算された欠陥モデルのディレクトリパス。
:param CAR_DIEL_el_dir: str 電子誘電率計算モデルのディレクトリパス、または電子誘電率の平均値。
:param CAR_DIEL_ion_dir: str イオン性誘電率計算モデルのディレクトリパス、またはイオン性誘電率の平均値。
:param CAR_charge_dir: str 欠陥の電荷を決定するための欠陥モデルのディレクトリパス。
:returns: None
"""
CAR_dir1 = CAR_IDIPOL_dir
CAR_dir2 = CAR_DIEL_el_dir
CAR_dir3 = CAR_DIEL_ion_dir
CAR_dir4 = CAR_charge_dir
vasp1 = tkVASP()
vasp2 = tkVASP()
vasp3 = tkVASP()
vasp4 = tkVASP()
base_path1 = vasp1.getdir(CAR_dir1)
base_path2 = vasp2.getdir(CAR_dir2)
base_path3 = vasp2.getdir(CAR_dir3)
base_path4 = vasp2.getdir(CAR_dir4)
print("")
print(f"Summary file: {summaryfile}")
INCAR_path1 = vasp1.get_INCAR(base_path1)
POSCAR_path1 = vasp1.get_POSCAR(base_path1)
POTCAR_path1 = vasp1.get_POTCAR(base_path1)
OUTCAR_path1 = vasp1.get_OUTCAR(base_path1)
print("Defect CAR dir(IDIPOL=4) : ", CAR_dir1)
print(" INCAR : ", INCAR_path1)
print(" POSCAR : ", POSCAR_path1)
print(" POTCAR : ", POTCAR_path1)
print(" OUTCAR : ", OUTCAR_path1)
INCAR_path2 = vasp2.get_INCAR(base_path2)
POSCAR_path2 = vasp2.get_POSCAR(base_path2)
OUTCAR_path2 = vasp2.get_OUTCAR(base_path2)
print("CAR dir(electron dielectric constant) : ", CAR_dir2)
print(" INCAR : ", INCAR_path2)
print(" POSCAR : ", POSCAR_path2)
print(" OUTCAR : ", OUTCAR_path2)
INCAR_path3 = vasp3.get_INCAR(base_path3)
POSCAR_path3 = vasp3.get_POSCAR(base_path3)
OUTCAR_path3 = vasp3.get_OUTCAR(base_path3)
print("CAR dir(ionic dielectric constant) : ", CAR_dir3)
print(" INCAR : ", INCAR_path3)
print(" POSCAR : ", POSCAR_path3)
print(" OUTCAR : ", OUTCAR_path3)
INCAR_path4 = vasp4.get_INCAR(base_path4)
POSCAR_path4 = vasp4.get_POSCAR(base_path4)
POTCAR_path4 = vasp4.get_POTCAR(base_path4)
OUTCAR_path4 = vasp4.get_OUTCAR(base_path4)
print("CAR dir(defect charge) : ", CAR_dir4)
print(" INCAR : ", INCAR_path4)
print(" POSCAR : ", POSCAR_path4)
print(" POTCAR : ", POTCAR_path4)
print(" OUTCAR : ", OUTCAR_path4)
print("")
print("##########################################")
print("# Defect model (IDIPOL=4)")
print("##########################################")
print("")
print("*** Read INCAR from [{}]".format(INCAR_path1))
incarinf1 = vasp1.read_incar_inf(INCAR_path1)
print("")
print("*** Read POTCAR from [{}]".format(POTCAR_path1))
potcarinf1 = vasp1.read_potcar_inf(POTCAR_path1)
nPPTypes = potcarinf1["nPPTypes"]
print("Pseudo potentials:")
print(" nPPTypes=", nPPTypes)
for i in range(nPPTypes):
ppname = potcarinf1["PP{}".format(i)]
atomname = potcarinf1["Atom{}".format(i)]
ne = potcarinf1["Ne{}".format(i)]
print(" {:2d}: {}".format(i, ppname))
print(" {} ne={}".format(atomname, ne))
"""
print("*** Read crystal structure from [{}]".format(POSCAR_path1))
cry1 = vasp1.read_poscar(POSCAR_path1)
if cry1 is None:
app.terminate(f"Error: Can not read [{POSCAR_path1}]", usage = usage, pause = True)
a1, b1, c1, alpha1, beta1, gamm1 = cry1.LatticeParameters()
Vcell = cry1.Volume()
cry1.PrintInf("cell")
print(" Vcell={} A^3".format(Vcell))
L = Vcell**(1.0/3.0)
print(" L={:10.6g} A".format(L))
poscarinf1 = vasp1.poscarinf
nAtomTypes = poscarinf1["nAtomTypes"]
print(" nAtomTypes=", nAtomTypes)
print(" AtomTypes =", poscarinf1["AtomTypes"])
print(" nAtomTypes=", poscarinf1["nAtomTypes"])
print(" nAtoms =", poscarinf1["nAtoms"])
"""
print("")
print("*** Read OUTCAR from [{}]".format(OUTCAR_path1))
outcarinf1 = vasp1.read_outcar_inf(OUTCAR_path1)
ISPIN = outcarinf1["ISPIN"]
NELECT = outcarinf1["NELECT"]
IDIPOL = outcarinf1.get("IDIPOL", None)
print("Information in OUTCAR:")
print(" ISPIN : ", ISPIN)
print(" IDIPOL: ", IDIPOL)
print(" NELECT: ", NELECT)
if IDIPOL is None or IDIPOL != 4:
app.terminate(f"Error: IDIPOL in {CAR_dir1} must be 4", usage = usage, pause = True)
print("")
print("##########################################")
print("# Ideal model (electronic dielectric constant)")
print("##########################################")
if IsFile(CAR_DIEL_el_dir) or IsDir(CAR_DIEL_el_dir):
print("*** Read electronic dielectric constants from [{}]".format(OUTCAR_path2))
dielinf_el = read_dielectric_inf_el(OUTCAR_path2)
print("")
print("*** Read INCAR from [{}]".format(INCAR_path2))
incarinf2 = vasp2.read_incar_inf(INCAR_path2)
print("*** Read OUTCAR from [{}]".format(OUTCAR_path2))
outcarinf2 = vasp3.read_outcar_inf(OUTCAR_path2)
ISPIN = outcarinf2["ISPIN"]
LOPTICS = outcarinf2.get("LOPTICS", None)
print("Information in OUTCAR:")
print(" ISPIN : ", ISPIN)
print(" LOPTICS: ", LOPTICS)
if LOPTICS is None or 't' not in LOPTICS.lower():
app.terminate(f"Error: LOPTICS in {CAR_dir2} must be .True.", usage = usage, pause = True)
eps_el = dielinf_el["eps_opt_r"]
print("e(el): |{:10.4g} {:10.4g} {:10.4g}|".format(eps_el[0][0], eps_el[0][1], eps_el[0][2]))
print(" |{:10.4g} {:10.4g} {:10.4g}|".format(eps_el[1][0], eps_el[1][1], eps_el[1][2]))
print(" |{:10.4g} {:10.4g} {:10.4g}|".format(eps_el[2][0], eps_el[2][1], eps_el[2][2]))
eps_av_el = (eps_el[0][0] + eps_el[1][1] + eps_el[2][2]) / 3.0
print(" trace average: {:10.6g}".format(eps_av_el))
else:
print("*** Read dielectric constant from argument [{}]".format(CAR_DIEL_el_dir))
eps_av_el = pfloat(CAR_DIEL_el_dir)
print(" e(el): {:10.6g}".format(eps_av_el))
print("")
print("##########################################")
print("# Ideal model (ionic dielectric constant)")
print("##########################################")
if IsFile(CAR_DIEL_ion_dir) or IsDir(CAR_DIEL_ion_dir):
print("*** Read ionic dielectric constants from [{}]".format(OUTCAR_path3))
dielinf = read_dielectric_inf_ion(OUTCAR_path3)
print("")
print("*** Read INCAR from [{}]".format(INCAR_path3))
incarinf3 = vasp3.read_incar_inf(INCAR_path3)
print("*** Read OUTCAR from [{}]".format(OUTCAR_path3))
outcarinf3 = vasp3.read_outcar_inf(OUTCAR_path3)
ISPIN = outcarinf3["ISPIN"]
LEPSILON = outcarinf3.get("LEPSILON", None)
print("Information in OUTCAR:")
print(" ISPIN : ", ISPIN)
print(" LEPSILON: ", LEPSILON)
if LEPSILON is None or 't' not in LEPSILON.lower():
app.terminate(f"Error: LEPSILON in {CAR_dir3} must be .True.", usage = usage, pause = True)
epstot = dielinf["epsilon(static/local field)"]
epsion = dielinf["epsilon(static/ion)"]
print("e(tot): |{:10.4g} {:10.4g} {:10.4g}|".format(epstot[0][0], epstot[0][1], epstot[0][2]))
print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epstot[1][0], epstot[1][1], epstot[1][2]))
print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epstot[2][0], epstot[2][1], epstot[2][2]))
epstotav = (epstot[0][0] + epstot[1][1] + epstot[2][2]) / 3.0
print(" trace average: {:10.6g}".format(epstotav))
print("e(ion): |{:10.4g} {:10.4g} {:10.4g}|".format(epsion[0][0], epsion[0][1], epsion[0][2]))
print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epsion[1][0], epsion[1][1], epsion[1][2]))
print(" |{:10.4g} {:10.4g} {:10.4g}|".format(epsion[2][0], epsion[2][1], epsion[2][2]))
epsionav = (epsion[0][0] + epsion[1][1] + epsion[2][2]) / 3.0
print(" trace average: {:10.6g}".format(epsionav))
else:
print("*** Read dielectric constant from argument [{}]".format(CAR_DIEL_ion_dir))
epsionav = pfloat(CAR_DIEL_ion_dir)
print(" eps(tot): {:10.6g}".format(epsionav))
eps_tot = epsionav + eps_av_el
print("")
print("*** Read IDIPOL information from [{}]".format(OUTCAR_path1))
dipoleinf = read_dipole_inf(OUTCAR_path1)
# print("d=", dipoleinf)
print(" dipolemoment:", dipoleinf["dipolemoment"])
print(" Tr[quadrupol] :", dipoleinf["Tr[quadrupol]"])
print(" Ecorr(di+quadrupol): {:12.6g} eV".format(dipoleinf["Ecorr(di+quadrupol)"]))
print(" Ecorr(with field) : {:12.6g} eV".format(dipoleinf["Ecorr(with field)"]))
print(" Ecorr(charge)=q^2 * aM / L : {:12.6g} eV".format(dipoleinf["Ecorr(charge)"]))
print("")
print("##########################################")
print("# Defect model (defect charge)")
print("##########################################")
print("")
print("*** Read crystal structure from [{}]".format(POSCAR_path4))
cry4 = vasp4.read_poscar(POSCAR_path4)
if cry4 is None:
app.terminate(f"Error: Can not read [{POSCAR_path4}]", usage = usage, pause = True)
a4, b4, c4, alpha4, beta4, gamm4 = cry4.LatticeParameters()
Vcell = cry4.Volume()
cry4.PrintInf("cell")
print(" Vcell={} A^3".format(Vcell))
L = Vcell**(1.0/3.0)
print(" L={:10.6g} A".format(L))
poscarinf4 = vasp4.poscarinf
nAtomTypes = poscarinf4["nAtomTypes"]
print(" nAtomTypes=", nAtomTypes)
print(" AtomTypes =", poscarinf4["AtomTypes"])
print(" nAtomTypes=", poscarinf4["nAtomTypes"])
print(" nAtoms =", poscarinf4["nAtoms"])
print("")
print("*** Read POTCAR from [{}]".format(POTCAR_path4))
potcarinf4 = vasp4.read_potcar_inf(POTCAR_path4)
nPPTypes = potcarinf4["nPPTypes"]
print("Pseudo potentials:")
print(" nPPTypes=", nPPTypes)
print("")
print("*** Read OUTCAR from [{}]".format(OUTCAR_path4))
outcarinf4 = vasp4.read_outcar_inf(OUTCAR_path4)
ISPIN = outcarinf4["ISPIN"]
NELECT = outcarinf4["NELECT"]
print("Information in OUTCAR:")
print(" ISPIN : ", ISPIN)
print(" NELECT: ", NELECT)
nVe = 0.0
for i in range(nAtomTypes):
key = "Ne{}".format(i)
ne = potcarinf4[key]
nAtoms = poscarinf4["nAtoms"][i]
nVe += nAtoms * ne
# print(" n=", i, ne, nAtoms, nVe)
print(" # of total valcence electrons: ", nVe)
q = nVe - NELECT
print(" NELECT =", NELECT)
print(" nVe =", nVe)
print(" q =", q)
print(" eps(tot)=", eps_tot, "e0")
print(" L =", L, " A")
print("")
print("Summary:")
print("Electronic dielectric constant (average): {:10.6g}".format(eps_av_el))
print("Ionic dielectric constant (average) : {:10.6g}".format(epsionav))
print("Total dielectric constant (average) : {:10.6g}".format(eps_tot))
if abs(q) <= 1.0e-3:
aM = dipoleinf["Ecorr(charge)"] * L
# print(" q^2 * aM : {:12.6g}".format(aM))
# dEcorr = 1.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / eps_tot
dEcorr = 0.0
print(f" q={q} is very small, regarded as zero")
print(" (1/3) * q^2 * aM / L / eps(total): {:12.6g} eV".format(dEcorr))
# dEcorr = 2.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / epstotav
# print(" (2/3) * q^2 * aM / L / eps(tot,av): {:12.6g} eV".format(dEcorr))
else:
aM = dipoleinf["Ecorr(charge)"] * L / q / q
print(" aM : {:12.6g} eV*A".format(aM))
dEcorr = 1.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / eps_tot
print(" (1/3) * q^2 * aM / L / eps(total): {:12.6g} eV".format(dEcorr))
# dEcorr = 2.0 / 3.0 * q * q * dipoleinf["Ecorr(charge)"] / epstotav
# print(" (2/3) * q^2 * aM / L / eps(tot,av): {:12.6g} eV".format(dEcorr))
print("")
print(f"Save summary to {summaryfile}")
ini = tkIniFile()
ini.write_from_scratch(summaryfile, 'results',
Car_dir1 = base_path1,
Car_dir2 = base_path2,
Car_dir3 = base_path3,
Car_dir4 = base_path4,
eps_av_el = eps_av_el,
epsionav = epsionav,
eps_tot = eps_tot,
q = q,
L = L,
Error_charge = dipoleinf["Ecorr(charge)"],
dEcorr = dEcorr
)
app.terminate("", usage = usage, pause = True)
[ドキュメント]
def main():
"""
スクリプトのメイン処理を実行する。
詳細説明:
この関数は、プログラムのエントリポイントとして機能します。
- グローバル変数をコマンドライン引数に基づいて更新します。
- ログファイルをリダイレクトして、標準出力とファイルに同時に記録します。
- 現在のモードに基づいて電荷補正計算関数 `correct_charged_defects` を呼び出します。
- 不正なモードが指定された場合はエラーメッセージを表示し、プログラムを終了します。
:returns: None
"""
global mode
global cifpath, poscarpath
global summaryfile
updatevars()
vasp = tkVASP()
base_path = vasp.getdir(CAR_IDIPOL_dir)
# logfile = app.replace_path(infile)
logfile = os.path.join(base_path, 'charge_correction-out.txt')
summaryfile = os.path.join(base_path, 'charge_correction-summary.prm')
print("")
print(f"Open logfile [{logfile}]")
app.redirect(targets = ["stdout", logfile], mode = 'w')
print("")
print("=============== Estimate the correction of charged defects and dipole interaction by VASP ============")
print("")
print("mode: ", mode)
if mode == 'Z':
correct_charged_defects(mode, CAR_IDIPOL_dir, CAR_DIEL_el_dir, CAR_DIEL_ion_dir, CAR_charge_dir)
else:
app.terminate(f"Error: Invalide mode [{mode}]", usage = usage, pause = True)
if __name__ == "__main__":
main()