draw_band_fs.py ダウンロード/コピー
draw_band_fs.py
draw_band_fs.py
1"""
2バンド構造、状態密度 (DOS)、フェルミエネルギー (EF)、有効質量 (m*)、およびフェルミ面を
3様々な次元とモデルでシミュレートしプロットするスクリプト。
4
5詳細説明:
6このスクリプトは、1次元から3次元までの単純な電子モデル（自由電子、タイトバインディング、
7またはほぼ自由電子）に基づき、以下の機能を提供する。
8- バンド構造の計算とプロット（1D直線パスまたは高対称点パスに沿って）。
9- DOSの計算とプロット（ガウスブロードニングまたはヒストグラム）。
10- 電子数 `nelec` からDOS積分を通じてフェルミエネルギー `EF` を決定する機能、
11 または `EF` を固定値として指定する機能。
12- 1Dバンドプロットに重ねて有効質量 `m*(k)` を表示する機能。
13- 1D、2D、3Dのフェルミ面をプロットする機能（3Dにはscikit-imageが必要）。
14
15利用可能なモデル:
16- 'free': 自由電子モデル。
17- 'tb': タイトバインディングモデル。
18- 'nfe': ほぼ自由電子 (NFE) モデル（逆格子ベクトル間のカップリングVを考慮）。
19
20関連リンク:
21:doc:`draw_band_fs_usage`
22"""
23#!/usr/bin/env python3
24# -*- coding: utf-8 -*-
25
26import numpy as np
27import matplotlib.pyplot as plt
28import matplotlib.colors as mcolors
29import argparse
30from mpl_toolkits.mplot3d.art3d import Poly3DCollection
31
32# skimage.measure is imported lazily in plot_3d_fermi_surface to avoid crash if not installed
33
34NORMALIZATION_UNIT = 2 * np.pi # k_int * 2π -> k_phys (a=1)
35
36
37# =============================================================================
38# Utilities
39# =============================================================================
40
41def reduce_to_first_bz(k_int: float) -> float:
42 """k_int を第1BZ [-0.5, 0.5) に還元（周期1）
43
44 詳細説明:
45 内部単位で表現されたk値を、周期境界条件を適用して第1ブリルアンゾーン
46 ([-0.5, 0.5) の範囲) にマッピングし直します。
47
48 :param k_int: 内部単位のk値。
49 :type k_int: float
50 :returns: 第1BZに還元されたk値。
51 :rtype: float
52 """
53 return ((k_int + 0.5) % 1.0) - 0.5
54
55
56def build_g_vecs(dim: int) -> np.ndarray:
57 """dim 次元の最小セットのG: 0 と ±e_i を物理単位で返す
58
59 詳細説明:
60 指定された次元における逆格子ベクトルGの最小セットを構築します。
61 これは、0ベクトルと、各軸方向の単位ベクトル (±e_i) を含みます。
62 結果は物理単位 (2π/a) に正規化されます。
63
64 :param dim: 空間次元 (1, 2, or 3)。
65 :type dim: int
66 :returns: Gベクトル群の配列。各行がGベクトル。
67 :rtype: np.ndarray
68 """
69 ints = [np.zeros(dim, dtype=float)]
70 for i in range(dim):
71 v = np.zeros(dim, dtype=float)
72 v[i] = 1.0
73 ints.append(v.copy())
74 ints.append(-v.copy())
75 return np.array(ints, dtype=float) * NORMALIZATION_UNIT
76
77
78def nband_of_model(dim: int, mode: str) -> int:
79 """このtoyモデルでのバンド本数（スピンは含めない）
80
81 詳細説明:
82 指定されたモデルと空間次元に基づき、スピン縮退を含まないバンドの総数を返します。
83 - 'tb': タイトバインディングモデルは1バンド。
84 - 'free'/'nfe': 自由電子モデルおよびほぼ自由電子モデルは、0ベクトルと±e_iベクトル
85 によって決定される数のバンド（1 + 2 * dim）を持ちます。
86
87 :param dim: 空間次元 (1, 2, or 3)。
88 :type dim: int
89 :param mode: モデルの種類 ('tb', 'free', 'nfe')。
90 :type mode: str
91 :returns: モデルのバンド本数。
92 :rtype: int
93 :raises ValueError: 未知のモデルが指定された場合。
94 """
95 if mode == "tb":
96 return 1
97 if mode in ("free", "nfe"):
98 return 1 + 2 * dim # G = 0, ±e_i
99 raise ValueError("bad mode")
100
101
102def determine_energy_range(E_samples: np.ndarray, E_range_user=None):
103 """エネルギー軸レンジを決める
104
105 詳細説明:
106 プロットやDOS/EF計算のためのエネルギー軸の最小値と最大値を決定します。
107 ユーザーが `E_range_user` を指定した場合はそれを採用し、
108 指定がない場合は、提供されたエネルギーサンプル `E_samples` の
109 最小値と最大値に5%のパディングを加えて自動的に決定します。
110
111 :param E_samples: サンプルされたエネルギー値の配列。
112 :type E_samples: np.ndarray
113 :param E_range_user: ユーザー指定のエネルギー範囲 (Emin, Emax)。Noneの場合は自動決定。
114 :type E_range_user: tuple[float, float], optional
115 :returns: 決定されたエネルギー範囲 (Emin, Emax)。
116 :rtype: tuple[float, float]
117 """
118 if E_range_user is not None:
119 Emin, Emax = map(float, E_range_user)
120 return Emin, Emax
121
122 Emin = float(np.min(E_samples))
123 Emax = float(np.max(E_samples))
124 pad = 0.05 * (Emax - Emin if Emax > Emin else 1.0)
125 return Emin - pad, Emax + pad
126
127
128# =============================================================================
129# Band energies at k
130# =============================================================================
131
132def get_band_energies_at_k_dim(k_vec_int: np.ndarray, model: str, dim: int, args) -> np.ndarray:
133 """指定されたk点におけるバンドエネルギーを計算する
134
135 詳細説明:
136 内部単位で与えられたkベクトル `k_vec_int` に対し、指定された `model`
137 と `dim` に基づいてバンドエネルギー（固有値）を計算します。
138 - 'tb': タイトバインディングモデル。
139 - 'free': 自由電子モデル。Gベクトルに対応する運動エネルギーを直接計算しソート。
140 - 'nfe': ほぼ自由電子モデル。Gベクトル間の周期ポテンシャル `args.v` を考慮した
141 ハミルトニアン行列を構築し、固有値を計算。
142
143 :param k_vec_int: 内部単位 (2π/a normalized) のkベクトル (dim,)。
144 :type k_vec_int: np.ndarray
145 :param model: モデルの種類 ('tb', 'free', 'nfe')。
146 :type model: str
147 :param dim: 空間次元 (1, 2, or 3)。
148 :type dim: int
149 :param args: コマンドライン引数オブジェクト (特に 'nfe' モデルの `args.v` に必要)。
150 :type args: argparse.Namespace
151 :returns: そのk点での固有エネルギーの配列 (Nband,)。
152 :rtype: np.ndarray
153 :raises ValueError: 未知のモデルが指定された場合。
154 """
155 k_phys = k_vec_int * NORMALIZATION_UNIT
156 G_VECS = build_g_vecs(dim)
157
158 if model == "tb":
159 return np.array([tb_energy(k_phys, dim)], dtype=float)
160
161 if model == "free":
162 # 自由電子モデル: 行列は作らず、各Gに対してエネルギーを直接計算
163 energies = []
164 for g in G_VECS:
165 dk = k_phys - g
166 # E = |k - G|^2 (定数係数は省略またはNORMALIZATION_UNITに内包)
167 energies.append(np.sum(dk**2))
168 return np.sort(np.array(energies, dtype=float))
169
170 if model == "nfe":
171 # ほぼ自由電子モデル: 逆格子ベクトル間のカップリングを考慮
172 num_G = len(G_VECS)
173 H = np.zeros((num_G, num_G), dtype=float)
174
175 # 対角成分: |k-G|^2
176 for i in range(num_G):
177 dk = k_phys - G_VECS[i]
178 H[i, i] = float(np.sum(dk**2))
179
180 # 非対角成分: 周期ポテンシャル V_G による散乱
181 for i in range(num_G):
182 for j in range(i + 1, num_G):
183 H[i, j] = args.v
184 H[j, i] = args.v
185
186 return np.linalg.eigvalsh(H)
187
188 raise ValueError(f"Unknown model: {model}")
189
190
191def tb_energy(k_phys: np.ndarray, dim: int) -> float:
192 """dim次元TB (t=1, onsite=0)
193
194 詳細説明:
195 指定された空間次元 `dim` におけるタイトバインディング (TB) モデルのエネルギーを計算します。
196 ホッピング積分 t=1、オンサイトエネルギーを0と仮定しています。
197 - 1D: -cos(k_x a)
198 - 2D: -(cos(k_x a) + cos(k_y a))
199 - 3D: -(cos(k_x a) + cos(k_y a) + cos(k_z a))
200 ここで k_phys は k*a の物理単位です。
201
202 :param k_phys: 物理単位のkベクトル。
203 :type k_phys: np.ndarray
204 :param dim: 空間次元 (1, 2, or 3)。
205 :type dim: int
206 :returns: 計算されたタイトバインディングエネルギー。
207 :rtype: float
208 :raises ValueError: `dim` が1, 2, 3以外の場合。
209 """
210 if dim == 1:
211 return -float(np.cos(k_phys[0]))
212 elif dim == 2:
213 return -float(np.cos(k_phys[0]) + np.cos(k_phys[1]))
214 elif dim == 3:
215 return -float(np.cos(k_phys[0]) + np.cos(k_phys[1]) + np.cos(k_phys[2]))
216 else:
217 raise ValueError("dim must be 1, 2, or 3")
218
219
220# =============================================================================
221# BZ sampling -> DOS -> EF(nelec)
222# =============================================================================
223
224def sample_energies_in_bz(args, dim: int) -> np.ndarray:
225 """dim 次元 BZ を格子サンプルして、全バンドのEを1次元配列で返す（状態列）
226
227 詳細説明:
228 指定された次元 `dim` のブリルアンゾーン内を `args.ef_res` で指定された
229 解像度で均一にサンプリングし、各k点における全バンドのエネルギーを計算します。
230 結果として得られるエネルギー値は1次元配列にフラット化されます。
231 `--mode tb` の場合は特別な高速パスを使用します。
232
233 :param args: コマンドライン引数オブジェクト (ef_res, modeなどに必要)。
234 :type args: argparse.Namespace
235 :param dim: 空間次元 (1, 2, or 3)。
236 :type dim: int
237 :returns: 計算された全エネルギー値の1次元配列。
238 :rtype: np.ndarray
239 :raises ValueError: `dim` が1, 2, 3以外の場合。
240 """
241 n = args.ef_res
242 k = np.linspace(-0.5, 0.5, n, endpoint=False)
243
244 if dim == 1:
245 pts = k[:, None]
246 elif dim == 2:
247 KX, KY = np.meshgrid(k, k, indexing="ij")
248 pts = np.stack([KX.ravel(), KY.ravel()], axis=1)
249 elif dim == 3:
250 KX, KY, KZ = np.meshgrid(k, k, k, indexing="ij")
251 pts = np.stack([KX.ravel(), KY.ravel(), KZ.ravel()], axis=1)
252 else:
253 raise ValueError("dim must be 1,2,3")
254
255 if args.mode == "tb":
256 k_phys = pts * NORMALIZATION_UNIT
257 if dim == 1:
258 E = -np.cos(k_phys[:, 0])
259 elif dim == 2:
260 E = -(np.cos(k_phys[:, 0]) + np.cos(k_phys[:, 1]))
261 elif dim == 3:
262 E = -(np.cos(k_phys[:, 0]) + np.cos(k_phys[:, 1]) + np.cos(k_phys[:, 2]))
263 else:
264 raise ValueError("dim must be 1,2,3")
265 return np.asarray(E, dtype=float)
266
267 evals_all = []
268 for kv in pts:
269 evals_all.append(get_band_energies_at_k_dim(kv, args.mode, dim, args))
270 return np.asarray(evals_all, dtype=float).ravel()
271
272
273def compute_dos(E_samples: np.ndarray, nE: int, sigma: float, E_range, scale_states: float):
274 """DOS をガウスブロードニングで作る
275
276 詳細説明:
277 与えられたエネルギーサンプル `E_samples` から、ガウスブロードニング
278 またはヒストグラム法を用いて状態密度 (DOS) を計算します。
279 `sigma` が0以下の場合はヒストグラムとして計算され、それ以外の場合は
280 ガウス関数を用いて平滑化されたDOSが計算されます。
281 結果のDOSは `scale_states` に正規化され、DOSの積分が `scale_states` に
282 ほぼ等しくなります。
283
284 :param E_samples: サンプルされたエネルギー値の配列。
285 :type E_samples: np.ndarray
286 :param nE: DOSを計算するためのエネルギーグリッドの点数。
287 :type nE: int
288 :param sigma: ガウスブロードニングの幅。0以下の場合はヒストグラムとして計算。
289 :type sigma: float
290 :param E_range: DOS計算を行うエネルギー範囲 (Emin, Emax)。
291 :type E_range: tuple[float, float]
292 :param scale_states: DOSをスケーリングする係数（通常はスピン縮退×バンド数）。
293 :type scale_states: float
294 :returns: エネルギーグリッドと対応するDOS値の配列。
295 :rtype: tuple[np.ndarray, np.ndarray]
296 """
297 E_samples = np.asarray(E_samples, dtype=float)
298 Emin, Emax = map(float, E_range)
299 E_grid = np.linspace(Emin, Emax, nE)
300
301 if sigma is None or sigma <= 0:
302 hist, edges = np.histogram(E_samples, bins=nE, range=(Emin, Emax), density=True)
303 centers = 0.5 * (edges[:-1] + edges[1:])
304 return centers, hist * scale_states
305
306 inv = 1.0 / (np.sqrt(2.0 * np.pi) * sigma)
307 dos = np.zeros_like(E_grid, dtype=float)
308
309 block = 20000
310 for i in range(0, len(E_samples), block):
311 Ei = E_samples[i:i + block]
312 x = (E_grid[:, None] - Ei[None, :]) / sigma
313 dos += inv * np.exp(-0.5 * x * x).sum(axis=1)
314
315 dos /= len(E_samples) # ∫dos dE ≈ 1
316 dos *= scale_states # ∫dos dE ≈ scale_states
317 return E_grid, dos
318
319
320def ef_from_dos(nelec: float, E_grid: np.ndarray, dos: np.ndarray, warn: bool = True) -> float:
321 """DOS を積分して N(E)=∫DOS dE を作り、N(EF)=nelec となる EF を内挿で求める。
322
323 詳細説明:
324 計算された状態密度 (DOS) をエネルギー `E_grid` に沿って積分し、
325 累積状態数 `N(E)` を構築します。
326 次に、指定された電子数 `nelec` に対応するフェルミエネルギー (EF) を
327 `N(E)` から線形内挿で求めます。
328 `nelec` がDOSの積分範囲外にある場合、EFは範囲の端にクランプされ、
329 `warn` がTrueの場合は警告メッセージが表示されます。
330
331 :param nelec: 単位胞あたりの電子数。
332 :type nelec: float
333 :param E_grid: DOSのエネルギーグリッドの配列。
334 :type E_grid: np.ndarray
335 :param dos: 対応するDOS値の配列。
336 :type dos: np.ndarray
337 :param warn: `nelec` がDOS範囲外の場合に警告を表示するかどうか。
338 :type warn: bool
339 :returns: 計算されたフェルミエネルギー (EF)。
340 :rtype: float
341 """
342 E_grid = np.asarray(E_grid, dtype=float)
343 dos = np.asarray(dos, dtype=float)
344
345 dE = np.diff(E_grid)
346 avg = 0.5 * (dos[:-1] + dos[1:])
347 Ncum = np.concatenate([[0.0], np.cumsum(avg * dE)])
348
349 # clamp target
350 nelec_clamped = float(np.clip(nelec, Ncum[0], Ncum[-1]))
351
352 if warn:
353 if nelec < Ncum[0] - 1e-12:
354 print(f"[WARN] nelec={nelec} is below integrated minimum ({Ncum[0]:.6g}). "
355 f"EF is clamped to Emin={E_grid[0]:.6g}.")
356 if nelec > Ncum[-1] + 1e-12:
357 print(f"[WARN] nelec={nelec} exceeds integrated states in DOS range ({Ncum[-1]:.6g}). "
358 f"EF is clamped to Emax={E_grid[-1]:.6g}. Consider widening --E_ef_range or omit it.")
359
360 EF = float(np.interp(nelec_clamped, Ncum, E_grid))
361 return EF
362
363
364# =============================================================================
365# kF along 1D path and root solving
366# =============================================================================
367
368def refine_root_bisection(func, a, b, tol=1e-10, max_iter=80):
369 """関数 func の根を二分法で求める
370
371 詳細説明:
372 与えられた関数 `func` の、区間 `[a, b]` 内での根を二分法アルゴリズムで探索します。
373 `func(a)` と `func(b)` の符号が異なる場合にのみ根が存在すると仮定します。
374 指定された許容誤差 `tol` または最大反復回数 `max_iter` に達するまで探索を続けます。
375
376 :param func: 根を探索する関数。`func(x)` は数値を返す必要があります。
377 :type func: callable
378 :param a: 探索区間の下限。
379 :type a: float
380 :param b: 探索区間の上限。
381 :type b: float
382 :param tol: 許容誤差。`abs(fmid)` または `abs(hi - lo)` がこの値以下になると探索を終了します。
383 :type tol: float, optional
384 :param max_iter: 最大反復回数。この回数を超えると探索を終了し、現在の区間の中間値を返します。
385 :type max_iter: int, optional
386 :returns: 根が見つかった場合はその値、`func(a)` と `func(b)` の符号が同じ場合はNone。
387 :rtype: float or None
388 """
389 fa = func(a)
390 fb = func(b)
391 if fa == 0.0:
392 return a
393 if fb == 0.0:
394 return b
395 if fa * fb > 0:
396 return None
397
398 lo, hi = a, b
399 flo = fa
400 for _ in range(max_iter):
401 mid = 0.5 * (lo + hi)
402 fmid = func(mid)
403 if abs(fmid) < tol or abs(hi - lo) < tol:
404 return mid
405 if flo * fmid <= 0:
406 hi = mid
407 else:
408 lo = mid
409 flo = fmid
410 return 0.5 * (lo + hi)
411
412
413def find_kf_crossings_1d(k_grid_int, energies_1d, EF, dim, args):
414 """1D直線パスでのkF点（直線プロット用）
415
416 詳細説明:
417 1Dパスに沿ったバンドエネルギー `energies_1d` をフェルミエネルギー `EF` と比較し、
418 `E(k) - EF = 0` となるkF点を二分法で探索します。
419 見つかったkF点は第1ブリルアンゾーンに還元され、`args.kf_merge_tol` で指定された
420 許容範囲内の近接するkF点はマージされます。
421
422 :param k_grid_int: 1Dパスのk点グリッド (内部単位)。
423 :type k_grid_int: np.ndarray
424 :param energies_1d: 各k点におけるバンドエネルギーの配列 (Nk, Nband)。
425 :type energies_1d: np.ndarray
426 :param EF: フェルミエネルギー。
427 :type EF: float
428 :param dim: 空間次元。バンドエネルギー計算のために使用されます。
429 :type dim: int
430 :param args: コマンドライン引数オブジェクト (kf_tol, kf_merge_tol, modeなどに必要)。
431 :type args: argparse.Namespace
432 :returns: (バンドインデックス, kF点) のタプルのリスト。
433 :rtype: list[tuple[int, float]]
434 """
435 Nk, Nb = energies_1d.shape
436 crossings = []
437
438 def get_Ek_band(kx_int, ib):
439 k_vec = np.zeros(dim, dtype=float)
440 k_vec[0] = kx_int
441 return float(get_band_energies_at_k_dim(k_vec, args.mode, dim, args)[ib])
442
443 for ib in range(Nb):
444 Ek = energies_1d[:, ib] - EF
445 for i in range(Nk - 1):
446 a = k_grid_int[i]
447 b = k_grid_int[i + 1]
448 fa = Ek[i]
449 fb = Ek[i + 1]
450 if fa == 0.0:
451 crossings.append((ib, a))
452 continue
453 if fa * fb > 0:
454 continue
455 root = refine_root_bisection(lambda x: get_Ek_band(x, ib) - EF, a, b, tol=args.kf_tol)
456 if root is not None:
457 crossings.append((ib, root))
458
459 # reduce to 1st BZ and merge per band
460 reduced = [(ib, reduce_to_first_bz(kx)) for (ib, kx) in crossings]
461 reduced.sort(key=lambda x: (x[0], x[1]))
462
463 merged = []
464 for ib, kx in reduced:
465 if not merged:
466 merged.append([ib, kx])
467 else:
468 ib0, k0 = merged[-1]
469 if ib == ib0 and abs(kx - k0) < args.kf_merge_tol:
470 merged[-1][1] = 0.5 * (k0 + kx)
471 else:
472 merged.append([ib, kx])
473
474 return [(ib, kx) for ib, kx in merged]
475
476
477# =============================================================================
478# Effective mass overlay
479# =============================================================================
480
481def effective_mass_1d(k_int: np.ndarray, E: np.ndarray) -> np.ndarray:
482 """1Dパス上で有効質量 m*(k) = 1 / (d^2E/d(ka)^2) を数値評価。
483
484 詳細説明:
485 1次元パスに沿ったバンドエネルギー `E` とk点 `k_int` を用いて、
486 電子の有効質量 `m*(k)` を数値的に計算します。
487 有効質量は運動量 `ka` に対するエネルギーの二階微分 `d^2E/d(ka)^2` の逆数として定義されます。
488 二階微分は中心差分法を用いて近似されます。k点数が少ない場合はNaNを返します。
489
490 :param k_int: 内部単位のk点配列。
491 :type k_int: np.ndarray
492 :param E: 各k点に対応するエネルギー配列。
493 :type E: np.ndarray
494 :returns: 計算された有効質量配列。計算できない点ではNaNが含まれます。
495 :rtype: np.ndarray
496 """
497 k_int = np.asarray(k_int, dtype=float)
498 E = np.asarray(E, dtype=float)
499
500 if len(k_int) < 3:
501 return np.full_like(E, np.nan, dtype=float)
502
503 dk = k_int[1] - k_int[0]
504 d2E = np.full_like(E, np.nan, dtype=float)
505 d2E[1:-1] = (E[2:] - 2.0 * E[1:-1] + E[:-2]) / (dk * dk)
506
507 with np.errstate(divide="ignore", invalid="ignore"):
508 mstar_kint = 1.0 / d2E
509
510 mstar_ka = mstar_kint * (NORMALIZATION_UNIT) ** 2 # 1/(d^2E/dk^2) * (dk/d(ka))^2 = 1/(d^2E/dk^2) * (1/ (2π/a)^2)^(-1) = 1/(d^2E/dk^2) * (2π/a)^2
511 return mstar_ka
512
513
514# =============================================================================
515# Band + DOS plotter (Straight path only)
516# =============================================================================
517
518def plot_band(args, dim: int, EF: float, E_samples: np.ndarray, E_range_plot, E_range_ef, Egrid_dos=None, dos=None):
519 """1Dバンドプロット、またはkx方向に沿った直線カット用 (DOS/m*対応)
520
521 詳細説明:
522 1次元のバンド構造、または高次元モデルのkx方向に沿った直線カットのバンド構造をプロットします。
523 フェルミエネルギー (EF) を破線で表示し、必要に応じてkF点を縦点線で示します。
524 `args.dos` がTrueの場合、DOSプロットをバンドプロットの隣に表示します。
525 `args.mstar` がTrueの場合、有効質量 `m*(k)` をバンドプロットの右軸に重ねて表示します。
526
527 :param args: コマンドライン引数オブジェクト。
528 :type args: argparse.Namespace
529 :param dim: 空間次元 (1, 2, or 3)。
530 :type dim: int
531 :param EF: フェルミエネルギー。
532 :type EF: float
533 :param E_samples: DOS計算に使用されたエネルギーサンプルの配列。DOSプロットがない場合はNoneでも可。
534 :type E_samples: np.ndarray or None
535 :param E_range_plot: プロットのY軸（エネルギー）範囲 (Emin, Emax)。
536 :type E_range_plot: tuple[float, float]
537 :param E_range_ef: DOS/EF計算に使用されたエネルギー範囲 (Emin, Emax)。DOSプロットがない場合はNoneでも可。
538 :type E_range_ef: tuple[float, float] or None
539 :param Egrid_dos: DOSのエネルギーグリッドの配列。`args.dos` がFalseの場合はNone。
540 :type Egrid_dos: np.ndarray or None
541 :param dos: 対応するDOS値の配列。`args.dos` がFalseの場合はNone。
542 :type dos: np.ndarray or None
543 :returns: なし。
544 :rtype: None
545 """
546 k_min, k_max = args.k_path_range
547 k_resolution = args.res * 3
548 k_range_int = np.linspace(k_min, k_max, k_resolution)
549
550 # 1D path energies
551 all_energies = []
552 for kx_int in k_range_int:
553 k_vec = np.zeros(dim, dtype=float)
554 k_vec[0] = kx_int
555 all_energies.append(get_band_energies_at_k_dim(k_vec, args.mode, dim, args))
556 all_energies = np.array(all_energies, dtype=float) # (Nk, Nband)
557
558 # kF roots on this cut
559 kf_list = find_kf_crossings_1d(k_range_int, all_energies, EF, dim, args)
560
561 Emin_plot, Emax_plot = E_range_plot
562
563 # DOS plot setup
564 is_dos_plot = args.dos and (Egrid_dos is not None)
565
566 if is_dos_plot:
567 Nband = nband_of_model(dim, args.mode)
568 scale_states = float(args.spin) * Nband
569 fig, (ax_band, ax_dos) = plt.subplots(1, 2, figsize=(11, 5), gridspec_kw={"wspace": 0.28})
570 else:
571 fig, ax_band = plt.subplots(figsize=(6.5, 5))
572 ax_dos = None
573
574 # --- band ---
575 is_reduced_zone = (k_min >= -0.5 and k_max <= 0.5)
576
577 for ib in range(all_energies.shape[1]):
578 ax_band.plot(k_range_int, all_energies[:, ib], linewidth=1.5, alpha=0.9)
579
580 if is_reduced_zone:
581 ax_band.set_xlim(-0.5, 0.5)
582
583 ax_band.set_ylim(Emin_plot, Emax_plot)
584 ax_band.set_title(rf"{args.type} ({args.mode}) $E_F={EF:.6g}$")
585 ax_band.axhline(y=EF, linestyle="--", linewidth=2, color="gray", label=r"$E_F$")
586
587 # kF vertical lines (reduced) + learned label
588 if kf_list:
589 for ib, kx_red in kf_list:
590 if k_min <= kx_red <= k_max:
591 ax_band.axvline(x=kx_red, linestyle=":", linewidth=1.2)
592 ax_band.plot([], [], linestyle=":", color="black", label=r"$k_F$")
593
594 # effective mass overlay
595 if args.mstar:
596 ax_m = ax_band.twinx()
597 for ib in range(all_energies.shape[1]):
598 mstar = effective_mass_1d(k_range_int, all_energies[:, ib])
599 mstar = np.clip(mstar, -5.0, 5.0) # required range
600 ax_m.plot(k_range_int, mstar, linestyle="--", linewidth=1.2, alpha=0.7, color="C1")
601
602 ax_m.set_ylabel(r"Effective mass $m^*(k)$")
603 ax_m.set_ylim(-5.0, 5.0)
604 ax_m.axhline(0.0, color="gray", linewidth=0.8, alpha=0.6)
605 ax_band.plot([], [], "C1--", label=r"$m^*(k)$")
606
607 ax_band.set_xlabel(r"$k_x / (2\pi/a)$")
608 ax_band.set_ylabel(r"Energy $E$")
609 ax_band.grid(True, linestyle="--", alpha=0.5)
610 ax_band.legend(loc="upper right")
611
612 # --- DOS ---
613 if ax_dos is not None:
614 ax_dos.plot(dos, Egrid_dos, linewidth=1.6)
615 ax_dos.axhline(EF, linestyle="--", linewidth=2, color="gray", label=r"$E_F$")
616
617 # DOSの縦軸はバンド表示レンジに合わせる（見た目の統一）
618 ax_dos.set_ylim(Emin_plot, Emax_plot)
619
620 ax_dos.set_xlabel(r"DOS (states / cell / energy)")
621 ax_dos.set_ylabel(r"Energy $E$")
622 ax_dos.set_title(rf"DOS ($\int$DOS dE = {scale_states:.3g})")
623 ax_dos.grid(True, linestyle="--", alpha=0.5)
624 ax_dos.legend(loc="upper right")
625
626 plt.tight_layout()
627 plt.show()
628
629
630# =============================================================================
631# High-symmetry K-Path Builder (for 2D/3D band plot)
632# =============================================================================
633
634def build_k_path(points, labels, res_per_segment=40):
635 """高対称点(points)をつなぐパスを生成する
636
637 詳細説明:
638 複数の高対称点 (`points`) を連続して結び、バンド図プロット用のKパスを構築します。
639 各パスセグメントは `res_per_segment` で指定された数のK点に分割されます。
640 K点ベクトル、Kパスに沿った累積距離、高対称点の位置とラベルを返します。
641
642 :param points: 高対称点の座標のリスト。例: `[[0.0, 0.0], [0.5, 0.0]]`。
643 :type points: list[list[float]]
644 :param labels: 各高対称点に対応するラベルのリスト。例: `["Gamma", "X"]`。
645 :type labels: list[str]
646 :param res_per_segment: 各パスセグメントあたりのK点解像度。
647 :type res_per_segment: int, optional
648 :returns:
649 - `k_vecs` (np.ndarray): パス上の全K点ベクトルの配列。
650 - `x_dist` (np.ndarray): パス上のK点に対応する累積距離の配列。
651 - `x_nodes` (list[float]): 高対称点に対応する累積距離のリスト。
652 - `labels` (list[str]): 高対称点ラベルのリスト。
653 :rtype: tuple[np.ndarray, np.ndarray, list[float], list[str]]
654 """
655 k_vecs_list = []
656 x_dist_list = []
657 x_nodes = []
658
659 current_dist = 0.0
660 x_nodes.append(current_dist)
661
662 # 修正済み: 始点を配列化してリストに追加
663 k_vecs_list.append(np.array(points[0]))
664 x_dist_list.append(np.array([current_dist]))
665
666 for i in range(len(points) - 1):
667 start = np.array(points[i])
668 end = np.array(points[i+1])
669
670 dist = np.linalg.norm(end - start)
671
672 # 点を生成（始点は重複するので含めず、終点を含める）
673 segment_k = np.linspace(start, end, res_per_segment + 1)[1:]
674 segment_dist = np.linspace(current_dist, current_dist + dist, res_per_segment + 1)[1:]
675
676 k_vecs_list.append(segment_k)
677 x_dist_list.append(segment_dist)
678
679 current_dist += dist
680 x_nodes.append(current_dist)
681
682 k_vecs = np.vstack(k_vecs_list)
683 x_dist = np.concatenate(x_dist_list)
684
685 return k_vecs, x_dist, x_nodes, labels
686
687
688def plot_band_with_path(args, dim: int, EF: float, E_range_plot):
689 """高対称点パスに沿ってバンド図を描画する (2D/3D兼用)
690
691 詳細説明:
692 2次元または3次元のバンド構造を高対称点パスに沿ってプロットします。
693 指定された次元に応じて、対応する標準的な高対称点（例: 2D正方格子のΓ-X-M-Γパス、
694 3D単純立方格子のΓ-X-M-Γ-R-Xパス）が定義され、そのパスに沿ってバンドエネルギーが計算されます。
695 フェルミエネルギー (EF) は破線で表示され、高対称点の位置は縦線とラベルで示されます。
696
697 :param args: コマンドライン引数オブジェクト。
698 :type args: argparse.Namespace
699 :param dim: 空間次元 (2 or 3)。
700 :type dim: int
701 :param EF: フェルミエネルギー。
702 :type EF: float
703 :param E_range_plot: プロットのY軸（エネルギー）範囲 (Emin, Emax)。
704 :type E_range_plot: tuple[float, float]
705 :returns: なし。
706 :rtype: None
707 """
708
709 # --- 1. 高対称点とパスの定義 ---
710 if dim == 2:
711 # 2D Square Lattice: Gamma -> X -> M -> Gamma
712 G_pt = [0.0, 0.0]
713 X_pt = [0.5, 0.0]
714 M_pt = [0.5, 0.5]
715 path_points = [G_pt, X_pt, M_pt, G_pt]
716 path_labels = [r"$\Gamma$", r"$X$", r"$M$", r"$\Gamma$"]
717
718 elif dim == 3:
719 # 3D Simple Cubic: Gamma -> X -> M -> Gamma -> R -> X
720 G_pt = [0.0, 0.0, 0.0]
721 X_pt = [0.5, 0.0, 0.0]
722 M_pt = [0.5, 0.5, 0.0]
723 R_pt = [0.5, 0.5, 0.5]
724 path_points = [G_pt, X_pt, M_pt, G_pt, R_pt, X_pt]
725 path_labels = [r"$\Gamma$", r"$X$", r"$M$", r"$\Gamma$", r"$R$", r"$X$"]
726
727 else:
728 # 1D band plotはplot_bandを使うため、通常ここは呼ばれない
729 print(f"[ERROR] Path plotting for dim={dim} is not supported.")
730 return
731
732 # --- 2. パス生成 ---
733 k_vecs, x_coords, node_pos, node_labs = build_k_path(path_points, path_labels, res_per_segment=args.res)
734
735 # --- 3. エネルギー計算 ---
736 all_energies = []
737 for kv in k_vecs:
738 es = get_band_energies_at_k_dim(kv, args.mode, dim, args)
739 all_energies.append(es)
740
741 all_energies = np.array(all_energies, dtype=float) # (Nk, Nband)
742
743 # --- 4. プロット ---
744 fig, ax = plt.subplots(figsize=(9, 6))
745
746 # バンド線
747 for ib in range(all_energies.shape[1]):
748 ax.plot(x_coords, all_energies[:, ib], linewidth=1.5, alpha=0.9, color="C0")
749
750 # フェルミ準位
751 ax.axhline(EF, linestyle="--", color="gray", linewidth=1.5, label=r"$E_F$")
752
753 # 高対称点の縦線とラベル
754 for x_pos in node_pos:
755 ax.axvline(x_pos, color="black", linestyle="-", linewidth=0.5)
756
757 ax.set_xticks(node_pos)
758 ax.set_xticklabels(node_labs, fontsize=12)
759
760 # グラフ設定
761 Emin, Emax = E_range_plot
762 ax.set_ylim(Emin, Emax)
763 ax.set_xlim(x_coords[0], x_coords[-1])
764 ax.set_ylabel("Energy")
765 ax.set_title(rf"{dim}D Band Structure ({args.mode}) $E_F={EF:.4g}$")
766 ax.grid(True, axis="y", linestyle="--", alpha=0.5)
767
768 plt.tight_layout()
769 plt.show()
770
771
772# =============================================================================
773# Fermi surface: build energy grid for FS plotting
774# =============================================================================
775
776def build_energy_grid_for_fs(args, dim: int):
777 """フェルミ面プロットのためにエネルギーグリッドを構築する
778
779 詳細説明:
780 指定された次元 `dim` において、ブリルアンゾーン全体を格子状にサンプリングし、
781 特定のバンドインデックス `args.fs_band` のエネルギー値を計算してグリッドを構築します。
782 モデル ('tb', 'free', 'nfe') に応じてエネルギー計算が実行されます。
783 このグリッドはフェルミ面の等値面または等高線プロットの基盤となります。
784
785 :param args: コマンドライン引数オブジェクト (res, fs_band, modeなどに必要)。
786 :type args: argparse.Namespace
787 :param dim: 空間次元 (1, 2, or 3)。
788 :type dim: int
789 :returns:
790 - `k_range_int` (np.ndarray): 各軸のk点範囲 (内部単位)。
791 - `E` (np.ndarray): 計算されたエネルギーグリッド。次元に応じた形状。
792 :rtype: tuple[np.ndarray, np.ndarray]
793 :raises ValueError: `dim` が1, 2, 3以外の場合。
794 """
795 N = args.res
796 k_range_int = np.linspace(-0.5, 0.5, N)
797 band_index = int(args.fs_band)
798
799 if dim == 1:
800 KX = k_range_int
801 E = np.zeros_like(KX, dtype=float)
802 for i, kx in enumerate(KX):
803 kv = np.array([kx], dtype=float)
804 evals = get_band_energies_at_k_dim(kv, args.mode, 1, args)
805 E[i] = evals[band_index]
806 return k_range_int, E
807
808 if dim == 2:
809 KX, KY = np.meshgrid(k_range_int, k_range_int, indexing="ij")
810 E = np.zeros_like(KX, dtype=float)
811
812 if args.mode == "tb":
813 E = -(np.cos(KX * NORMALIZATION_UNIT) + np.cos(KY * NORMALIZATION_UNIT))
814 return k_range_int, E
815
816 if args.mode == "free":
817 G = build_g_vecs(2)
818 KX_phys = KX * NORMALIZATION_UNIT
819 KY_phys = KY * NORMALIZATION_UNIT
820 dists = []
821 for g in G:
822 d2 = (KX_phys - g[0]) ** 2 + (KY_phys - g[1]) ** 2
823 dists.append(d2)
824 all_bands = np.stack(dists, axis=0)
825 all_bands.sort(axis=0)
826 if band_index < len(G): E = all_bands[band_index]
827 else: E = all_bands[-1]
828 return k_range_int, E
829
830 if args.mode == "nfe":
831 it = np.nditer([KX, KY], flags=["multi_index"])
832 for kx, ky in it:
833 kv = np.array([float(kx), float(ky)], dtype=float)
834 evals = get_band_energies_at_k_dim(kv, "nfe", 2, args)
835 E[it.multi_index] = evals[band_index]
836 return k_range_int, E
837
838 if dim == 3:
839 KX, KY, KZ = np.meshgrid(k_range_int, k_range_int, k_range_int, indexing="ij")
840 E = np.zeros_like(KX, dtype=float)
841
842 if args.mode == "tb":
843 E = -(np.cos(KX * NORMALIZATION_UNIT) + np.cos(KY * NORMALIZATION_UNIT) + np.cos(KZ * NORMALIZATION_UNIT))
844 return k_range_int, E
845
846 if args.mode == "free":
847 G = build_g_vecs(3)
848 KX_phys = KX * NORMALIZATION_UNIT
849 KY_phys = KY * NORMALIZATION_UNIT
850 KZ_phys = KZ * NORMALIZATION_UNIT
851 dists = []
852 for g in G:
853 d2 = (KX_phys - g[0]) ** 2 + (KY_phys - g[1]) ** 2 + (KZ_phys - g[2]) ** 2
854 dists.append(d2)
855 all_bands = np.stack(dists, axis=0)
856 all_bands.sort(axis=0)
857 if band_index < len(G): E = all_bands[band_index]
858 else: E = all_bands[-1]
859 return k_range_int, E
860
861 if args.mode == "nfe":
862 it = np.nditer([KX, KY, KZ], flags=["multi_index"])
863 for kx, ky, kz in it:
864 kv = np.array([float(kx), float(ky), float(kz)], dtype=float)
865 evals = get_band_energies_at_k_dim(kv, "nfe", 3, args)
866 E[it.multi_index] = evals[band_index]
867 return k_range_int, E
868
869 raise ValueError("dim must be 1,2,3")
870
871
872# =============================================================================
873# Fermi surface plotters (1dfs/2dfs/3dfs)
874# =============================================================================
875
876def plot_1d_fermi_points(k_range_int: np.ndarray, E_line: np.ndarray, EF: float, args):
877 """1次元のフェルミ点（kF）をプロットする
878
879 詳細説明:
880 1次元のバンドエネルギー `E_line` とフェルミエネルギー `EF` を用いて、
881 `E(k) = EF` となるkF点をプロットします。
882 kF点はバンドとフェルミ準位の交点として見つけられ、第1ブリルアンゾーンに還元されて表示されます。
883
884 :param k_range_int: k点範囲の配列 (内部単位)。
885 :type k_range_int: np.ndarray
886 :param E_line: 各k点に対応するバンドエネルギーの配列。
887 :type E_line: np.ndarray
888 :param EF: フェルミエネルギー。
889 :type EF: float
890 :param args: コマンドライン引数オブジェクト (modeなどに必要)。
891 :type args: argparse.Namespace
892 :returns: なし。
893 :rtype: None
894 """
895 fig, ax = plt.subplots(figsize=(7, 4))
896 ax.plot(k_range_int, E_line, linewidth=2.0)
897 ax.axhline(EF, linestyle="--", linewidth=2, label=r"$E_F$")
898 roots = []
899 Ek = E_line - EF
900 for i in range(len(k_range_int) - 1):
901 a, b = k_range_int[i], k_range_int[i + 1]
902 fa, fb = Ek[i], Ek[i + 1]
903 if fa == 0: roots.append(a)
904 if fa * fb > 0: continue
905 if fb != fa:
906 r = a + (b - a) * (-fa) / (fb - fa)
907 roots.append(r)
908 roots_red = sorted({round(reduce_to_first_bz(r), 12) for r in roots})
909 for r in roots_red: ax.axvline(r, linestyle=":", linewidth=1.5)
910 if roots_red: ax.plot([], [], linestyle=":", label=r"$k_F$")
911 ax.set_xlabel(r"$k / (2\pi/a)$")
912 ax.set_ylabel(r"Energy $E$")
913 ax.set_title(rf"1D Fermi points ({args.mode}) $E_F={EF:.6g}$")
914 ax.grid(True, linestyle="--", alpha=0.5)
915 ax.legend(loc="best")
916 plt.tight_layout()
917 plt.show()
918
919
920def plot_2d_fermi_surface(E_grid_2d: np.ndarray, EF: float, args, k_range_int: np.ndarray):
921 """2次元フェルミ面をコンタープロットとして表示する
922
923 詳細説明:
924 2次元のエネルギーグリッド `E_grid_2d` を用いて、フェルミエネルギー `EF` の
925 等高線を2Dのフェルミ面としてプロットします。
926 `E <= EF` の領域は占有領域として塗りつぶされ、
927 `args.nx_view` に基づいて複数ブリルアンゾーンにわたって表示されるため、
928 周期的な性質が視覚化されます。
929
930 :param E_grid_2d: 2次元エネルギーグリッドの配列。
931 :type E_grid_2d: np.ndarray
932 :param EF: フェルミエネルギー。
933 :type EF: float
934 :param args: コマンドライン引数オブジェクト (nx_view, modeなどに必要)。
935 :type args: argparse.Namespace
936 :param k_range_int: k点範囲の配列 (内部単位)。グリッドの軸生成に使用。
937 :type k_range_int: np.ndarray
938 :returns: なし。
939 :rtype: None
940 """
941 view_limit = args.nx_view
942 N_TILES = int(np.ceil(2 * view_limit / 0.5))
943 if N_TILES % 2 == 0: N_TILES += 1
944 E_tile = np.tile(E_grid_2d, (N_TILES, N_TILES))
945 k_min_total = -0.5 * N_TILES
946 k_max_total = 0.5 * N_TILES
947 k_range_ext = np.linspace(k_min_total, k_max_total, len(k_range_int) * N_TILES)
948 fig, ax = plt.subplots(figsize=(8, 8))
949 E_min = float(np.min(E_tile))
950 ax.contourf(k_range_ext, k_range_ext, E_tile.T, levels=[E_min, EF], colors=["lightblue"], alpha=0.6, extend="neither")
951 ax.contour(k_range_ext, k_range_ext, E_tile.T, levels=[EF], colors=["red"], linewidths=2)
952 limit_int = 0.5
953 ax.set_xlim(-view_limit, view_limit)
954 ax.set_ylim(-view_limit, view_limit)
955 ax.plot([limit_int, limit_int], [-view_limit, view_limit], "k--", lw=1)
956 ax.plot([-limit_int, -limit_int], [-view_limit, view_limit], "k--", lw=1)
957 ax.plot([-view_limit, view_limit], [limit_int, limit_int], "k--", lw=1)
958 ax.plot([-view_limit, view_limit], [-limit_int, -limit_int], "k--", lw=1)
959 ax.set_aspect("equal")
960 ax.set_xlabel(r"$k_x / (2\pi/a)$")
961 ax.set_ylabel(r"$k_y / (2\pi/a)$")
962 ax.set_title(rf"2D Fermi Surface (Periodic, {args.mode}) $E_F={EF:.6g}$")
963 ax.grid(True, linestyle="--", alpha=0.5)
964 ax.fill_between([], [], color="lightblue", alpha=0.6, label=r"Occupied ($E \leq E_F$)")
965 ax.legend(loc="lower left")
966 plt.tight_layout()
967 plt.show()
968
969
970def plot_3d_fermi_surface(E_grid_3d: np.ndarray, EF: float, args, k_range_int: np.ndarray):
971 """3次元フェルミ面を3Dサーフェスとしてプロットする
972
973 詳細説明:
974 3次元のエネルギーグリッド `E_grid_3d` を用いて、フェルミエネルギー `EF` の
975 等値面（フェルミ面）をMarching Cubesアルゴリズムで抽出し、3Dグラフィックとして表示します。
976 scikit-imageライブラリが必要です。面には光源を考慮したシェーディングが適用され、
977 `args.fs_alpha` で透明度を調整できます。
978
979 :param E_grid_3d: 3次元エネルギーグリッドの配列。
980 :type E_grid_3d: np.ndarray
981 :param EF: フェルミエネルギー。
982 :type EF: float
983 :param args: コマンドライン引数オブジェクト (mode, fs_alphaなどに必要)。
984 :type args: argparse.Namespace
985 :param k_range_int: k点範囲の配列 (内部単位)。グリッドの座標変換に使用。
986 :type k_range_int: np.ndarray
987 :returns: なし。
988 :rtype: None
989 """
990 print("[INFO] Marching Cubes...")
991 try:
992 from skimage import measure
993 except ImportError:
994 print("\n[ERROR] 'scikit-image' is required for 3D Fermi surface plotting.")
995 print("Please install it via: pip install scikit-image\n")
996 return
997
998 try:
999 verts, faces, normals, values = measure.marching_cubes(E_grid_3d, EF, spacing=(1, 1, 1))
1000 step_size = 1.0 / (len(k_range_int) - 1)
1001 verts = verts * step_size - 0.5 # [0..N] -> [-0.5..0.5]
1002 except ValueError:
1003 print("[ERROR] Isosurface could not be generated at this EF.")
1004 return
1005
1006 fig = plt.figure(figsize=(10, 8))
1007 ax = fig.add_subplot(111, projection="3d")
1008 ax.set_box_aspect([1, 1, 1])
1009
1010 tris = verts[faces]
1011 v1 = tris[:, 1] - tris[:, 0]
1012 v2 = tris[:, 2] - tris[:, 0]
1013 face_normals = np.cross(v1, v2)
1014 norm_mag = np.linalg.norm(face_normals, axis=1)
1015 norm_mag[norm_mag == 0] = 1.0
1016 face_normals = face_normals / norm_mag[:, np.newaxis]
1017
1018 light_A = np.array([1.0, 1.0, 1.0])
1019 light_A /= np.linalg.norm(light_A)
1020 base_color_A = "cyan" if args.mode != "tb" else "orange"
1021 rgb_A = np.array(mcolors.to_rgb(base_color_A))
1022 base_color_B = "navy" if args.mode != "tb" else "darkred"
1023 rgb_B = np.array(mcolors.to_rgb(base_color_B))
1024
1025 dots = np.dot(face_normals, light_A)
1026 weights = (dots + 1.0) / 2.0
1027 weights = np.clip(weights, 0.0, 1.0)[:, np.newaxis]
1028
1029 face_colors_rgb = weights * rgb_A + (1.0 - weights) * rgb_B
1030 alpha_val = args.fs_alpha
1031 face_colors_rgba = np.hstack((face_colors_rgb, np.full((len(faces), 1), alpha_val)))
1032
1033 mesh = Poly3DCollection(verts[faces], facecolors=face_colors_rgba)
1034 mesh.set_edgecolor("black")
1035 mesh.set_linewidth(0.05)
1036 ax.add_collection3d(mesh)
1037
1038 limit = 0.5
1039 ax.set_xlim(-limit, limit)
1040 ax.set_ylim(-limit, limit)
1041 ax.set_zlim(-limit, limit)
1042
1043 title_str = rf"3D Fermi Surface ({args.mode}) $E_F={EF:.6g}$" + "\n" + r"(Dual Light Shading)"
1044 ax.set_title(title_str)
1045 ax.set_xlabel(r"$k_x / (2\pi/a)$")
1046 ax.set_ylabel(r"$k_y / (2\pi/a)$")
1047 ax.set_zlabel(r"$k_z / (2\pi/a)$")
1048 plt.show()
1049
1050
1051# =============================================================================
1052# Dispatcher (main function)
1053# =============================================================================
1054
1055def dim_from_type(t: str) -> int:
1056 """プロットタイプ文字列から空間次元を決定する
1057
1058 詳細説明:
1059 コマンドライン引数で与えられたプロットタイプ文字列に基づいて、
1060 対応する空間次元（1, 2, または 3）を返します。
1061 例: "1dband" -> 1, "2dfs" -> 2。
1062
1063 :param t: プロットタイプを示す文字列 ('1dband', '2dband', '3dband', '1dfs', '2dfs', '3dfs')。
1064 :type t: str
1065 :returns: 対応する空間次元。
1066 :rtype: int
1067 :raises ValueError: 未知のプロットタイプが指定された場合。
1068 """
1069 if t in ("1dband", "1dfs"): return 1
1070 if t in ("2dband", "2dfs"): return 2
1071 if t in ("3dband", "3dfs"): return 3
1072 raise ValueError("bad --type")
1073
1074
1075def main():
1076 """スクリプトのメインエントリポイント
1077
1078 詳細説明:
1079 コマンドライン引数をパースし、指定された `type` (バンド図、DOS、フェルミ面) に基づいて
1080 適切な計算とプロットを実行します。
1081 フェルミエネルギー (EF) は、電子数 `nelec` と計算されたDOSから決定されるか、
1082 `--ef_fixed` オプションで固定値を指定することもできます。
1083 各プロットタイプ (`1dband`, `2dband`, `3dband`, `1dfs`, `2dfs`, `3dfs`) に応じて、
1084 対応するプロット関数を呼び出します。
1085 """
1086 parser = argparse.ArgumentParser(
1087 description="Toy band / DOS / EF(nelec from DOS integral) / kF / m*(k) and FS (1dfs/2dfs/3dfs) simulator."
1088 )
1089
1090 parser.add_argument("--type", type=str, default="3dband",
1091 choices=["1dband", "2dband", "3dband", "1dfs", "2dfs", "3dfs"],
1092 help="Output type: band plots or Fermi surface plots.")
1093
1094 parser.add_argument("--mode", type=str, default="tb",
1095 choices=["free", "tb", "nfe"],
1096 help="Model: free / tb / nfe (toy NFE Hamiltonian).")
1097
1098 # --- 新規追加オプション ---
1099 parser.add_argument("--ef_fixed", type=float, default=None,
1100 help="Fixed Fermi energy (EF). If set, skips DOS calculation (nelec is ignored).")
1101 # --------------------------
1102
1103 # electron count
1104 parser.add_argument("--nelec", type=float, default=1.0,
1105 help="Electrons per unit cell (including spin). EF is obtained from DOS integral (ignored if --ef_fixed is set).")
1106 parser.add_argument("--spin", type=float, default=2.0,
1107 help="Spin degeneracy factor (default 2).")
1108
1109 # grid / resolution
1110 parser.add_argument("--res", type=int, default=40,
1111 help="Resolution for FS grids / band path segments.")
1112 parser.add_argument("--ef_res", type=int, default=24,
1113 help="BZ sampling resolution per axis for DOS/EF (cost grows as ef_res^dim; ignored if --ef_fixed is set).")
1114
1115 # NFE coupling
1116 parser.add_argument("--v", type=float, default=2.0,
1117 help="Off-diagonal coupling strength in toy NFE Hamiltonian.")
1118
1119 # band path
1120 parser.add_argument("--k_path_range", type=float, nargs=2, default=[-0.5, 0.5],
1121 metavar=("KMIN", "KMAX"),
1122 help="kx range for 1D band/1D cut plots (internal unit).")
1123
1124 # plot Y-range (visual only)
1125 parser.add_argument("--E_range", type=float, nargs=2, default=None, metavar=("EMIN", "EMAX"),
1126 help="Energy axis range for plotting (visual only).")
1127
1128 # EF/DOS calculation range (important!)
1129 parser.add_argument("--E_ef_range", type=float, nargs=2, default=None, metavar=("EMIN", "EMAX"),
1130 help="Energy range used for DOS/EF calculation. Default: auto from sampled energies. (ignored if --ef_fixed is set).")
1131
1132 # DOS options
1133 parser.add_argument("--dos", action="store_true",
1134 help="Also plot DOS (next to band). Requires DOS calculation (--ef_fixed must be None).")
1135 parser.add_argument("--dos_nE", type=int, default=800,
1136 help="Energy grid points for DOS.")
1137 parser.add_argument("--dos_sigma", type=float, default=0.15,
1138 help="Gaussian broadening sigma for DOS (<=0 => histogram).")
1139
1140 # kF roots controls
1141 parser.add_argument("--kf_tol", type=float, default=1e-10,
1142 help="Tolerance for kF root finding.")
1143 parser.add_argument("--kf_merge_tol", type=float, default=1e-4,
1144 help="Merge tolerance for reduced kF lines (internal unit).")
1145
1146 # effective mass overlay
1147 parser.add_argument("--mstar", action="store_true",
1148 help="Overlay effective mass m*(k) on band plot (right axis), clipped to [-5,5].")
1149
1150 # FS options
1151 parser.add_argument("--nx_view", type=float, default=1.5,
1152 help="2D FS periodic view range (internal unit).")
1153 parser.add_argument("--fs_band", type=int, default=0,
1154 help="Band index used for FS isosurface/contour (default 0).")
1155 parser.add_argument("--fs_alpha", type=float, default=0.6,
1156 help="Alpha for 3D FS surface.")
1157
1158 args = parser.parse_args()
1159 dim = dim_from_type(args.type)
1160
1161 # basic checks
1162 if args.spin <= 0: raise ValueError("--spin must be > 0")
1163 if args.nelec < 0: raise ValueError("--nelec must be >= 0")
1164
1165 # --- EF/DOS計算ロジック分岐 ---
1166 E_samples = None
1167 Egrid_dos = None
1168 dos = None
1169
1170 if args.ef_fixed is not None:
1171 # EFが固定値の場合
1172 EF = float(args.ef_fixed)
1173 print(f"[INFO] EF (Fixed) = {EF:.6g}")
1174 if args.dos:
1175 print("[WARN] --dos ignored because --ef_fixed is set (skipping time-consuming BZ sampling).")
1176 args.dos = False # DOSプロットを強制的にオフ
1177
1178 # プロットレンジ決定のために、最小限のエネルギー範囲を計算する（ここでは省略し、プロットレンジはユーザー指定かデフォルトに依存）
1179 if args.E_range is None:
1180 print("[WARN] EF is fixed but --E_range is not set. Using default auto-range logic might be inaccurate.")
1181
1182 E_range_ef = None # DOS計算はしないので不要
1183
1184 else:
1185 # nelecからEFを計算する場合
1186 print(f"[INFO] Sampling energies in {dim}D BZ for DOS/EF (ef_res={args.ef_res}) ...")
1187 E_samples = sample_energies_in_bz(args, dim)
1188
1189 Nband = nband_of_model(dim, args.mode)
1190 scale_states = float(args.spin) * Nband
1191 if args.nelec > scale_states + 1e-12:
1192 print(f"[WARN] nelec={args.nelec} exceeds max capacity spin*Nband={scale_states:.6g} "
1193 f"for this toy basis; EF will saturate near top of sampled states.")
1194
1195 E_range_ef = determine_energy_range(E_samples, args.E_ef_range)
1196
1197 # DOS and EF from DOS integral
1198 Egrid_dos, dos = compute_dos(
1199 E_samples,
1200 nE=args.dos_nE,
1201 sigma=args.dos_sigma,
1202 E_range=E_range_ef,
1203 scale_states=scale_states
1204 )
1205 EF = ef_from_dos(args.nelec, Egrid_dos, dos, warn=True)
1206 print(f"[INFO] EF(from DOS integral) = {EF:.6g} (nelec={args.nelec}, spin={args.spin}, Nband={Nband})")
1207
1208 # ------------------------------------
1209
1210 # プロット用レンジ決定（固定値/自動計算に関わらず）
1211 if args.E_range is not None:
1212 E_range_plot = (float(args.E_range[0]), float(args.E_range[1]))
1213 elif E_range_ef is not None:
1214 E_range_plot = E_range_ef
1215 else:
1216 # EF固定時にE_range指定がない場合のフォールバック（DOS計算はしないので、EF周辺を適当に表示）
1217 pad = 2.0
1218 E_range_plot = (EF - pad, EF + pad)
1219
1220 # 5) dispatch
1221 if args.type in ("2dband", "3dband"):
1222 plot_band_with_path(args, dim, EF, E_range_plot)
1223 return 0
1224
1225 if args.type == "1dband":
1226 # 1Dは直線プロット（DOS情報が必要なら渡す）
1227 plot_band(args, dim, EF, E_samples, E_range_plot, E_range_ef, Egrid_dos, dos)
1228 return 0
1229
1230 # FS energy grids and plots use EF computed above
1231 if args.type == "1dfs":
1232 k_range_int, E_line = build_energy_grid_for_fs(args, 1)
1233 plot_1d_fermi_points(k_range_int, E_line, EF, args)
1234 return 0
1235
1236 if args.type == "2dfs":
1237 k_range_int, E2 = build_energy_grid_for_fs(args, 2)
1238 plot_2d_fermi_surface(E2, EF, args, k_range_int)
1239 return 0
1240
1241 if args.type == "3dfs":
1242 k_range_int, E3 = build_energy_grid_for_fs(args, 3)
1243 plot_3d_fermi_surface(E3, EF, args, k_range_int)
1244 return 0
1245
1246 raise ValueError("unknown --type")
1247
1248
1249if __name__ == "__main__":
1250 raise SystemExit(main())