vasp_correction_bandfilling.py ダウンロード/コピー
vasp_correction_bandfilling.py をダウンロード
vasp_correction_bandfilling.py
vasp_correction_bandfilling.py
1"""
2概要:
3 VASP計算結果から欠陥モデルのバンドフィリング補正を推定します。
4
5詳細説明:
6 本スクリプトは、欠陥モデルの形成エネルギー計算において必要となるバンドフィリング補正 dE_BF を算出します。
7 理想結晶と欠陥結晶のVASP計算結果 (OUTCAR, DOSCAR, EIGENVAL) を比較し、
8 dEVBMによるバンドアラインメント後の状態占有率から電子と正孔によるエネルギー補正量を評価します。
9 また、計算結果のDOSとバンド構造をグラフとして可視化し、サマリーファイルに結果を保存します。
10
11関連リンク:
12 vasp_correction_bandfilling_usage
13"""
14import os
15import sys
16import shutil
17import glob
18import csv
19import re
20import numpy as np
21from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, pi
22from scipy.interpolate import interp1d
23from pprint import pprint
24from matplotlib import pyplot as plt
25
26
27from tklib.tkfile import tkFile
28from tklib.tkinifile import tkIniFile
29import tklib.tkre as tkre
30from tklib.tkutils import save_csv
31from tklib.tkapplication import tkApplication
32from tklib.tkutils import IsDir, IsFile, SplitFilePath
33from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
34from tklib.tksci.tksci import Reduce01, Round
35from tklib.tksci.tkconvolution import convolution, convolve_func
36from tklib.tksci.tkconvolution import convolve_xydata
37from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
38from tklib.tkcrystal.tkcif import tkCIF, tkCIFData
39from tklib.tkcrystal.tkcrystal import tkCrystal
40from tklib.tkcrystal.tkatomtype import tkAtomType
41from tklib.tkcrystal.tkvasp import tkVASP
42from tklib.tkcrystal.tkbandstructure import plot_band
43from tklib.tksci import tkequation
44import tklib.tkcsv
45
46
47#================================
48# global parameters
49#================================
50debug = 0
51
52validmodes = ['BF']
53mode = 'BF'
54
55CAR_dir1 = '.'
56CAR_dir2 = '.'
57summaryfile = None
58
59# dEVBM obtained by VBM correction
60dEVBM = ''
61
62# FWHM of Gaussian smearing function for DOS
63WG_DOS = 0.1 # eV
64
65
66# Energy to search edge energies
67EF0 = 0.1
68# occupancy threthold to separate HOMO and LUMO
69occ_th = 0.5
70
71# Critera DOS value to search band edges
72DOSth = 1.0e-5
73
74# Max number of electrons occupies one state
75Nemax = 1.0
76
77# Plot configuration
78band_marker_size = 4
79band_marker_edge_width = 0.5
80
81Emin = -10.0 # eV
82Emax = 10.0 # eV
83
84#colors = ['#000000', '#ff0000', '#00aa00', '#0000ff', '#aaaa00', '#ff00ff', '#00ffff', '#aa0000', '#00aa00', '#0000aa']
85#colors = ['k', 'r', 'g', 'b', 'y', 'm', 'c']
86figsize = (10, 8)
87fontsize = 16
88titlefontsize = 12
89labelfontsize = 12
90legend_fontsize = 12
91
92
93app = tkApplication()
94
95
96#=============================
97# Treat argments
98#=============================
99def usage():
100 """
101 概要:
102 スクリプトの正しい使用方法を標準出力に表示します。
103
104 詳細説明:
105 スクリプトの実行に必要なコマンドライン引数のフォーマットと、
106 各引数の意味について説明します。
107
108 戻り値:
109 :returns: なし
110 :rtype: None
111 """
112 global mode
113
114 if mode not in validmodes:
115 mode = validmodes[0]
116
117 print("")
118 print("Usage:")
119 print(" (a) python {} mode CAR_dir(ideal) CAR_dir(defect) dEVBM WG Emin Emax EF0".format(sys.argv[0]))
120 print(" mode: {}".format(validmodes))
121 print(" ex: python {} {} {} {} {} {} {} {} {}".format(sys.argv[0], mode, CAR_dir1, CAR_dir2, dEVBM, WG_DOS, Emin, Emax, EF0))
122
123def updatevars():
124 """
125 概要:
126 コマンドライン引数に基づいてグローバル変数を更新します。
127
128 詳細説明:
129 sys.argvから引数を読み込み、mode、CAR_dir1、CAR_dir2、dEVBM、WG_DOS、Emin、Emax、EF0
130 などのグローバルパラメータに値を設定します。不正なモードが指定された場合はエラーで終了します。
131
132 戻り値:
133 :returns: なし
134 :rtype: None
135 """
136 global mode
137 global CAR_dir1, CAR_dir2
138 global EF0, dEVBM
139 global Emin, Emax, WG_DOS
140
141 mode = getarg (1, mode)
142 CAR_dir1 = getarg (2, CAR_dir1)
143 CAR_dir2 = getarg (3, CAR_dir2)
144 dEVBM = getarg (4, dEVBM)
145 WG_DOS = getfloatarg(5, WG_DOS)
146 Emin = getfloatarg(6, Emin)
147 Emax = getfloatarg(7, Emax)
148 EF0 = getfloatarg(8, EF0)
149 if mode == 'BF':
150 pass
151 else:
152 terminate("Error: Invalide mode [{}]".format(mode), usage = usage)
153
154
155#def save_csv(path, headerlist, datalist, is_print = 0):
156
157def BF_correction(mode, CAR_path1, CAR_path2):
158 """
159 概要:
160 バンドフィリング補正を計算し、結果をファイルに保存し、グラフをプロットします。
161
162 詳細説明:
163 理想結晶と欠陥結晶のVASP計算結果ファイル（POSCAR, OUTCAR, DOSCAR, EIGENVAL）を読み込み、
164 それぞれのバンド構造、DOS、フェルミ準位などの情報を取得します。
165 dEVBMを用いて欠陥結晶のエネルギー準位を補正した後、バンドフィリング補正量（dEh, dEe, dEtot）を算出します。
166 計算結果はサマリーファイルに保存され、DOSとバンド構造のグラフが生成されます。
167
168 引数:
169 :param mode: 実行モード。現在 BF のみがサポートされています。
170 :type mode: str
171 :param CAR_path1: 理想結晶のVASP計算結果ディレクトリへのパス。
172 :type CAR_path1: str
173 :param CAR_path2: 欠陥結晶のVASP計算結果ディレクトリへのパス。
174 :type CAR_path2: str
175
176 戻り値:
177 :returns: なし
178 :rtype: None
179 """
180 global dEVBM, EF0
181
182 debug = 0
183
184 vasp1 = tkVASP()
185 vasp2 = tkVASP()
186
187 base_path1 = vasp1.getdir(CAR_path1)
188 base_path2 = vasp2.getdir(CAR_path2)
189
190 print("")
191 print(f"Summary file: {summaryfile}")
192 INCAR_path1 = vasp1.get_INCAR(base_path1)
193 POSCAR_path1 = vasp1.get_POSCAR(base_path1)
194 OUTCAR_path1 = vasp1.get_OUTCAR(base_path1)
195 DOSCAR_path1 = vasp1.get_VASPPath(base_path1, 'DOSCAR')
196 EIGENVAL_path1 = vasp1.get_VASPPath(base_path1, 'EIGENVAL')
197 INCAR_path2 = vasp2.get_INCAR(base_path2)
198 POSCAR_path2 = vasp2.get_POSCAR(base_path2)
199 OUTCAR_path2 = vasp2.get_OUTCAR(base_path2)
200 DOSCAR_path2 = vasp2.get_VASPPath(base_path2, 'DOSCAR')
201 EIGENVAL_path2 = vasp2.get_VASPPath(base_path2, 'EIGENVAL')
202 print("CAR dir(ideal) : ", CAR_path1)
203 print(" INCAR : ", INCAR_path1)
204 print(" POSCAR : ", POSCAR_path1)
205 print(" OUTCAR : ", OUTCAR_path1)
206 print(" DOSCAR : ", DOSCAR_path1)
207 print(" EIGENVAL: ", EIGENVAL_path1)
208 print("CAR dir(defect): ", CAR_path2)
209 print(" INCAR : ", INCAR_path2)
210 print(" POSCAR : ", POSCAR_path2)
211 print(" OUTCAR : ", OUTCAR_path2)
212 print(" DOSCAR : ", DOSCAR_path2)
213 print(" EIGENVAL: ", EIGENVAL_path2)
214
215 print("")
216 print("dEVBM: {} eV".format(dEVBM))
217 print("")
218 print("DOS plot")
219 print(f" FWHM of Gauss function for convolution: {WG_DOS} eV")
220# print("")
221# print("EF0: {} eV".format(EF0))
222
223 dEVBM = pfloat(dEVBM, '')
224 if dEVBM == '':
225 app.terminate("Error: dEVBM value must be given", pause = True)
226
227 print("")
228 print("Ideal crystal model:")
229 print("Read crystal structure from [{}]".format(POSCAR_path1))
230 cry1 = vasp1.read_poscar(POSCAR_path1)
231 if cry1 is None:
232 terminate("Error: Can not read [{}]".format(POSCAR_path1), usage = usage)
233
234 a1, b1, c1, alpha1, beta1, gamm1 = cry1.LatticeParameters()
235 Vcell1 = cry1.Volume()
236 types1 = cry1.AtomTypeList()
237 ntypes1 = len(types1)
238 sites1 = cry1.ExpandedAtomSiteList()
239 nsites1 = len(sites1)
240 atomnamelist1 = cry1.get_atom_name_list(mode = 'all', NameOnly = True)
241 poslist1 = cry1.get_position_list(mode = 'all', IsReduce01 = True)
242 cry1.PrintInf("cell")
243
244 print("")
245 print("Read OUTCAR from [{}]".format(OUTCAR_path1))
246 outcarinf1 = vasp1.read_outcar_inf(OUTCAR_path1)
247 if outcarinf1 is None:
248 terminate("Error: Can not read [{}]".format(OUTCAR_path1), usage = usage)
249 EF1 = vasp1.outcarinf["EF"]
250 ISPIN1 = outcarinf1["ISPIN"]
251# ETOT1 = outcarinf1["TOTEN"]
252 ETOT1 = outcarinf1["TOTEN_zero_sigma"]
253 print(f"ISPIN={ISPIN1}")
254 print(f"ETOT(sigma->0)={ETOT1}")
255
256 print("")
257 print("Defect crystal model:")
258 print("Read from [{}]".format(POSCAR_path2))
259 cry2 = vasp2.read_poscar(POSCAR_path2)
260 if cry2 is None:
261 terminate("Error: Can not read [{}]".format(POSCAR_path2), usage = usage)
262
263 a2, b2, c2, alpha2, beta2, gamm2 = cry2.LatticeParameters()
264 Vcell2 = cry2.Volume()
265 types2 = cry2.AtomTypeList()
266 ntypes2 = len(types2)
267 sites2 = cry2.ExpandedAtomSiteList()
268 nsites2 = len(sites2)
269 atomnamelist2 = cry2.get_atom_name_list(mode = 'all', NameOnly = True)
270 poslist2 = cry2.get_position_list(mode = 'all', IsReduce01 = True)
271 cry2.PrintInf("cell")
272
273 print("")
274 print("Read OUTCAR from [{}]".format(OUTCAR_path2))
275 outcarinf2 = vasp2.read_outcar_inf(OUTCAR_path2)
276 if outcarinf2 is None:
277 terminate("Error: Can not read [{}]".format(OUTCAR_path2), usage = usage)
278 EF2 = vasp2.outcarinf["EF"]
279 ISPIN2 = outcarinf2["ISPIN"]
280# ETOT2 = outcarinf2["TOTEN"]
281 ETOT2 = outcarinf2["TOTEN_zero_sigma"]
282 print(f"ISPIN={ISPIN2}")
283 print(f"ETOT(sigma->0)={ETOT2}")
284
285# Estimate supercell size
286 print("")
287 nx = int(a2 / a1 + 0.2)
288 ny = int(b2 / b1 + 0.2)
289 nz = int(c2 / c1 + 0.2)
290 print("Supercell multipliers of the defect crystal with respect to the ideal crystal")
291 print(" (nx, ny, nz) = ({}, {}, {})".format(nx, ny, nz))
292
293 print("")
294 print(f"Ideal crystal model: Read DOSCAR from [{DOSCAR_path1}]")
295 doscarinf1 = vasp1.read_doscar(DOSCAR_path1, EF = EF1, normalize_E = False, unit = '/cm3')
296 nE1 = doscarinf1["nE"]
297 E1 = doscarinf1["E"]
298 if nE1 == 0:
299 Edmin1 = 0.0
300 Edmax1 = 0.0
301 else:
302 Edmin1 = min(E1)
303 Edmax1 = max(E1)
304 tDOSup1 = doscarinf1["TotalDOSup"]
305 Neup1 = doscarinf1["Neup"]
306 tDOSdn1 = doscarinf1["TotalDOSdn"]
307 Nedn1 = doscarinf1["Nedn"]
308 if nE1 > 0:
309 print(f"EF={EF1} eV")
310 print(f"DOS E range: {Edmin1:10.6g} - {Edmax1:10.6g} eV, {nE1} points")
311 print(f" Energy is not normalized")
312 else:
313 print("")
314 app.terminate(f"Error: Zero nE for DOS(E) in {DOSCAR_path1}", pause = True)
315
316 print(f"Read EIGENVAL from [{DOSCAR_path1}]")
317 eigenvalinf1 = vasp1.read_eigenval(EIGENVAL_path1, EF = EF1, normalize_E = False)
318 nk1 = eigenvalinf1["nk"]
319 nLevels1 = eigenvalinf1["nLevels"]
320 bandedgeinf1a = vasp1.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf1, ISPIN = ISPIN1, occ_th = occ_th)
321 bandedgeinf1b = vasp1.gbandedges(base_path1, OUTCAR_path1, EIGENVAL_path1, EF1)
322 print("k points in EIGENVAL:")
323 print("nk=", nk1)
324 print("nLevels=", nLevels1)
325 """
326 for i in range(nk1):
327 el = eigenvalinf1['EList'][i]
328 kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el
329 print(" ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g} {:8.4f} {:12.4f}".format(kx, ky, kz, wk, dk, ktot))
330 """
331
332 EV1a = bandedgeinf1a["EV"]
333 EC1a = bandedgeinf1a["EC"]
334 Eg1a = bandedgeinf1a["Eg"]
335 EV1 = bandedgeinf1b["EVBM"]
336 EC1 = bandedgeinf1b["ECBM"]
337 Eg1 = bandedgeinf1b["Eg"]
338 EHOMO1 = bandedgeinf1a["EHOMO"]
339 ELUMO1 = bandedgeinf1a["ELUMO"]
340 print(f"Band edge from EIGENVAL:")
341 print(f"EF0={EF0:10.6f} eV")
342 print(f" find_band_edges: EV={EV1a:10.6f} EC={EC1a:10.6f} Eg={Eg1a:10.6f} eV")
343 print(f" HOMO:{EHOMO1:10.6f} eV")
344 print(f" LUMO:{ELUMO1:10.6f} eV")
345 print(f" gbandedges() : EV={EV1:10.6f} EC={EC1:10.6f} Eg={Eg1:10.6f} eV")
346
347 print("")
348 print(f"Defect crystal model: Read DOSCAR from [{DOSCAR_path2}]")
349 doscarinf2 = vasp2.read_doscar(DOSCAR_path2, EF = EF2, normalize_E = False, unit = '/cm3')
350 nE2 = doscarinf2["nE"]
351 E2 = doscarinf2["E"]
352 if nE2 == 0:
353 Edmin2 = 0.0
354 Edmax2 = 0.0
355 else:
356 Edmin2 = min(E2)
357 Edmax2 = max(E2)
358 tDOSup2 = doscarinf2["TotalDOSup"]
359 Neup2 = doscarinf2["Neup"]
360 tDOSdn2 = doscarinf2["TotalDOSdn"]
361 Nedn2 = doscarinf2["Nedn"]
362 if nE2 > 0:
363 print(f"EF={EF2} eV")
364 print(f"DOS E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points")
365 print(f" Energy is not normalized")
366
367 print(f"Read EIGENVAL from [{DOSCAR_path2}]")
368 eigenvalinf2 = vasp2.read_eigenval(EIGENVAL_path2, EF = EF2, normalize_E = False)
369 nk2 = eigenvalinf2["nk"]
370 nLevels2 = eigenvalinf2["nLevels"]
371 bandedgeinf2a = vasp2.find_band_edges_from_eigenval(EF0 = EF0, eigenvalinf = eigenvalinf2, ISPIN = ISPIN2, occ_th = occ_th)
372 bandedgeinf2b = vasp2.gbandedges(base_path2, OUTCAR_path2, EIGENVAL_path2, EF2)
373 print("k points in EIGENVAL:")
374 print("nk=", nk2)
375 print("nLevels=", nLevels2)
376 for i in range(nk2):
377 el = eigenvalinf2['EList'][i]
378 kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = el
379 print(" ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g} {:8.4f} {:12.4f}".format(kx, ky, kz, wk, dk, ktot))
380
381 EV2a = bandedgeinf1a["EV"]
382 EC2a = bandedgeinf1a["EC"]
383 Eg2a = bandedgeinf1a["Eg"]
384 EV2 = bandedgeinf1b["EVBM"]
385 EC2 = bandedgeinf1b["ECBM"]
386 Eg2 = bandedgeinf1b["Eg"]
387 EHOMO2 = bandedgeinf1a["EHOMO"]
388 ELUMO2 = bandedgeinf1a["ELUMO"]
389 print(f"Band edge from EIGENVAL:")
390 print(f"EF0={EF0:10.6f} eV")
391 print(f" find_band_edges: EV={EV2a:10.6f} EC={EC2a:10.6f} Eg={Eg2a:10.6f} eV")
392 print(f" HOMO:{EHOMO2:10.6f} eV")
393 print(f" LUMO:{ELUMO2:10.6f} eV")
394 print(f" gbandedges() : EV={EV2:10.6f} EC={EC2:10.6f} Eg={Eg2:10.6f} eV")
395
396 print("")
397 print("DOS plot range")
398 print(f" Ideal crystal model:")
399 print(f" EF={EF1} eV")
400 print(f" E range: {Edmin1:10.6g} - {Edmax1:10.6g} eV, {nE1} points")
401 """
402 print(f" DOS(E) is scaled for the size of the defect model by {nx}x{ny}x{nz}")
403 k = nx * ny * nz * 1.0
404 for i in range(len(E1)):
405 tDOSup1[i] *= k
406 Neup1[i] *= k
407 if ISPIN1 == 2:
408 tDOSdn2[i] *= k
409 Nedn1[i] *= k
410 """
411
412 print(f" Defect crystal model:")
413 print(f" Original EF={EF2} eV")
414 print(f" Original E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points")
415 EF2 -= dEVBM
416 Edmin2 -= dEVBM
417 Edmax2 -= dEVBM
418 for i in range(len(E2)):
419 E2[i] -= dEVBM
420 for ik in range(nk2):
421 el = eigenvalinf2['EList'][ik]
422 kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el]
423 Eups[ik] -= dEVBM
424 if ISPIN2 == 2:
425 Edns[ik] -= dEVBM
426 print(f" Adjusted EF={EF2} eV")
427 print(f" Adjusted E range: {Edmin2:10.6g} - {Edmax2:10.6g} eV, {nE2} points")
428
429 print("")
430 print("Calculate band filling energy")
431 totalw = 0.0
432 totalNe = 0.0
433 dEe = 0.0
434 dEh = 0.0
435 for ik in range(nk2):
436 el = eigenvalinf2['EList'][ik]
437 kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el]
438 totalw += wk
439 print(" k[{:3d}]: ({:8.4f} {:8.4f} {:8.4f}) w={:8.4g}".format(ik, kx, ky, kz, wk))
440 for iE in range(nLevels2):
441 Eup = Eups[iE]
442 neup = occups[iE]
443 totalNe += neup * wk
444
445 if Eup <= EV1 and abs(neup - Nemax) > 1.0e-6:
446 dEh += wk * (Nemax - neup) * (Eup - EV1)
447 print(f" dEh: up iE={iE:3d} E={Eup:10.6g} eV dNh={Nemax - neup:10.4g} dEh={dEh:12.6g} eV")
448 if Eup >= EC1 and neup > 1.0e-6:
449 dEe += wk * neup * (Eup - EC1)
450 print(f" dEe: up iE={iE:3d} E={Eup:10.6g} eV dNe={neup:10.4g} dEe={dEe:12.6g} eV")
451 if ISPIN2 == 2:
452 Edn = Edns[iE]
453 nedn = occdns[iE]
454 totalNe += nedn * wk
455
456 if Edn <= EV1 and abs(nedn - Nemax) > 1.0e-6:
457 dEh += wk * (Nemax - nedn) * (Edn - EV1)
458 print(f" dEh: dn iE={iE:3d} E={Edn:10.6g} eV dNh={Nemax - nedn:10.4g} dEh={dEh:12.6g} eV")
459 if Edn >= EC1 and nedn > 1.0e-6:
460 dEe += wk * nedn * (Edn - EC1)
461 print(f" dEe: dn iE={iE:3d} E={Edn:10.6g} eV dNe={nedn:10.4g} dEe={dEe:12.6g} eV")
462
463 print("")
464 print("From defect crystal model:")
465 print(" Total weight: ", totalw)
466 print(" NELECT : ", outcarinf2["NELECT"])
467 dEtot = dEh + dEe
468 if ISPIN2 == 1:
469 totalNe *= 2.0
470 dEh *= 2.0
471 dEe *= 2.0
472 dEtot *= 2.0
473 print(" Count spin degeneracy of 2:")
474 print(" Total Ne : ", totalNe)
475 print(" dEh = {:12.6g} eV".format(dEh))
476 print(" dEe = {:12.6g} eV".format(dEe))
477 print(" dEtot = {:12.6g} eV".format(dEtot))
478
479 print("")
480 print(f"Save summary to {summaryfile}")
481 ini = tkIniFile()
482 ini.write_from_scratch(summaryfile, 'results',
483 Car_dir1 = base_path1,
484 Car_dir2 = base_path2,
485 Etot1 = ETOT1 * nx * ny * nz,
486 Etot2 = ETOT2,
487 totalNe = totalNe,
488 dEh = dEh,
489 dEe = dEe,
490 dEtot = dEtot
491 )
492
493
494#=============================
495# Plot graphs
496#=============================
497 print("")
498
499 fig = plt.figure(figsize = figsize)
500 ax1 = fig.add_subplot(2, 2, 1)
501 ax2 = fig.add_subplot(2, 2, 3)
502 axband = fig.add_subplot(1, 2, 2)
503
504 print("E1=", E1)
505 tDOSup_s1 = convolution(E1, tDOSup1, WG_DOS, func_type = 'gauss')
506 ax1.plot(E1, tDOSup_s1, label = 'up(ideal)', linestyle = '-', linewidth = 0.5, color = 'black')
507 if ISPIN1 == 2:
508 tDOSdn_s1 = convolution(E, tDOSdn1, WG_DOS, func_type = 'gauss')
509 ax1.plot(E1, tDOSdn_s1, label = 'dn(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'black')
510 ylim = ax1.get_ylim()
511 ax1.plot([EV1, EV1], ylim, label = '$E_V$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red')
512 ax1.plot([EC1, EC1], ylim, label = '$E_C$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red')
513 ax1.plot([EF1, EF1], ylim, label = '$E_F$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'blue')
514# ax1.plot([EF0, EF0], ylim, label = '$E_{F0}$', linestyle = 'dashed', linewidth = 0.5, color = 'blue')
515# ax1.plot([EHOMO1, EHOMO1], ylim, label = '$E_{HOMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
516# ax1.plot([ELUMO1, ELUMO1], ylim, label = '$E_{LUMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
517# ax1.set_xlabel("$E$ (eV)", fontsize = fontsize)
518 ax1.set_ylabel("DOS (states/cm$^3$/eV)", fontsize = fontsize)
519# ax1.set_xlim([Emin, Emax])
520 ax1.tick_params(labelsize = fontsize)
521
522 tDOSup_s2 = convolution(E2, tDOSup2, WG_DOS, func_type = 'gauss')
523 ax1.plot(E2, tDOSup_s2, label = 'up(defect)', linestyle = '-', linewidth = 1.0, color = 'red')
524 if ISPIN2 == 2:
525 tDOSdn_s2 = convolution(E, tDOSdn2, WG_DOS, func_type = 'gauss')
526 ax1.plot(E2, tDOSdn_s2, label = 'dn(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'red')
527 ax1.plot([EF2, EF2], ylim, label = '$E_F$(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'blue')
528 ax1.tick_params(labelsize = fontsize)
529 ax1.set_title(f"{CAR_path1}\n{CAR_path2}", fontsize = titlefontsize)
530 ax1.legend(fontsize = legend_fontsize)
531
532 ax2.plot(E2, Neup2, label = 'up(defect)', linestyle = '-', linewidth = 0.5, color = 'black')
533 if ISPIN2 == 2:
534 ax2.plot(E2, Nedn2, label = 'dn(defect)', linestyle = '-', linewidth = 0.5, color = 'red')
535 ylim = ax2.get_ylim()
536 ax2.plot([EV1, EV1], ylim, label = '$E_V$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red')
537 ax2.plot([EC1, EC1], ylim, label = '$E_C$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'red')
538# ax2.plot([EF0, EF0], ylim, label = '$E_{F0}$', linestyle = 'dashed', linewidth = 0.5, color = 'blue')
539# ax2.plot([EHOMO1, EHOMO1], ylim, label = '$E_{HOMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
540# ax2.plot([ELUMO1, ELUMO1], ylim, label = '$E_{LUMO}$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'green')
541 ax2.plot([EF1, EF1], ylim, label = '$E_F$(ideal)', linestyle = 'dashed', linewidth = 0.5, color = 'blue')
542 ax2.plot([EF2, EF2], ylim, label = '$E_F$(defect)', linestyle = 'dashed', linewidth = 1.0, color = 'blue')
543 ax2.set_xlabel("$E$ (eV)", fontsize = fontsize)
544 ax2.set_ylabel("$N_e$ (states/unit cell)", fontsize = fontsize)
545# ax2.set_xlim([Emin, Emax])
546 ax2.legend(fontsize = legend_fontsize)
547 ax2.tick_params(labelsize = fontsize)
548
549 xk = []
550 yEup = []
551 yNup = []
552 yEdn = []
553 yNdn = []
554 for i in range(nk2):
555 el = eigenvalinf2['EList'][i]
556 kx, ky, kz, wk, dk, ktot, Eups, occups, Edns, occdns = [*el]
557 xk.append(ktot)
558 yEup.append(Eups)
559 yNup.append(occups)
560 if ISPIN2 == 2:
561 yEdn.append(Edns)
562 yNdn.append(occdns)
563
564 plot_band(plt, axband, ISPIN2, xk, yEup, yEdn, [Emin, Emax], occups = yNup, occdns = yNdn,
565 ktotallist = None, ktotal_namelist = None,
566 EV = EV1, EC = EC1, EF = EF2, EHOMO = None, ELUMO = None,
567 EFlabel = '$E_{F}$(defect)',
568 markersize = band_marker_size, markeredgewidth = band_marker_edge_width)
569 axband.set_title(f"dEh={dEh:12.6g} eV\ndEe={dEe:12.6g} eV\ndEtot={dEtot:12.6g} eV", fontsize = titlefontsize)
570
571 plt.tight_layout()
572
573
574# Rearange the graph axes so that they are not overlapped
575 plt.tight_layout()
576
577 """
578 plt.show()
579 print("")
580 print("Close graph window to terminate")
581 """
582
583 plt.pause(0.1)
584
585 app.terminate("", pause = True)
586
587
588def main():
589 """
590 概要:
591 スクリプトのメインエントリポイント。引数の処理とバンドフィリング補正の実行を行います。
592
593 詳細説明:
594 コマンドライン引数をパースしてグローバル変数を更新し、ログファイルをセットアップします。
595 指定されたモードに応じてBF_correction関数を呼び出し、バンドフィリング補正計算を実行します。
596
597 戻り値:
598 :returns: なし
599 :rtype: None
600 """
601 global mode
602 global cifpath, poscarpath
603 global summaryfile
604
605 updatevars()
606
607 vasp = tkVASP()
608 base_path = vasp.getdir(CAR_dir2)
609 logfile = os.path.join(base_path, 'BF_correction-out.txt')
610 summaryfile = os.path.join(base_path, 'BF_correction-summary.prm')
611 print("")
612 print(f"Open logfile [{logfile}]")
613 app.redirect(targets = ["stdout", logfile], mode = 'w')
614
615 print("")
616 print("=============== Estimate band filling correction for defect model calculated by VASP ============")
617 print("")
618 print("mode: ", mode)
619
620 if mode == 'BF':
621 BF_correction(mode, CAR_dir1, CAR_dir2)
622 else:
623 terminate(f"Error: Invalide mode [{mode}]", usage = usage)
624
625
626if __name__ == "__main__":
627 main()