bz_draw.py ダウンロード/コピー
bz_draw.py
bz_draw.py
1#!/usr/bin/env python3
2"""
3第一ブリルアンゾーンまたは実空間のウィグナー・ザイツセル描画スクリプト。
4
5このスクリプトは、指定された結晶学的パラメータと格子タイプに基づいて、
6第一ブリルアンゾーン（逆空間）またはウィグナー・ザイツセル（実空間）を3Dで可視化します。
7さらに、周囲の周期的なセルをタイリング表示する機能、高対称点のパス表示（逆空間の場合）、
8およびインタラクティブなパラメータ調整機能（スライダーとApplyボタン）を提供します。
9
10:doc:`bz_draw_usage`
11"""
12import argparse
13import numpy as np
14import sys, math, functools
15import matplotlib.pyplot as plt
16from matplotlib.widgets import Slider, Button
17from mpl_toolkits.mplot3d.art3d import Poly3DCollection
18from scipy.spatial import Voronoi
19import warnings
20warnings.filterwarnings(
21 "ignore",
22 category=DeprecationWarning,
23 module=r"spglib\.spglib",
24 message=r"dict interface \(SpglibDataset\['.*'\]\) is deprecated"
25)
26
27# ================== 幾何ユーティリティ ==================
28def build_conventional_basis(a, b, c, alpha, beta, gamma):
29 """
30 従来の結晶学的パラメータから直交座標系における格子基底ベクトルを構築します。
31
32 最初の基底ベクトルをx軸に沿わせ、2番目の基底ベクトルをxy平面に置き、
33 3番目の基底ベクトルを正のz成分を持つように配置します。
34
35 :param a: 格子定数 a (Å)
36 :type a: float
37 :param b: 格子定数 b (Å)
38 :type b: float
39 :param c: 格子定数 c (Å)
40 :type c: float
41 :param alpha: 角度 alpha (度)
42 :type alpha: float
43 :param beta: 角度 beta (度)
44 :type beta: float
45 :param gamma: 角度 gamma (度)
46 :type gamma: float
47 :returns: 構築された3x3の格子基底行列。各列が基底ベクトル (a_vec, b_vec, c_vec) を表します。
48 :rtype: numpy.ndarray
49 """
50 ar, br, gr = np.radians([alpha, beta, gamma])
51 a_vec = np.array([a, 0.0, 0.0])
52 b_vec = np.array([b*np.cos(gr), b*np.sin(gr), 0.0])
53 cx = c*np.cos(br)
54 denom = np.sin(gr) if abs(np.sin(gr)) > 1e-12 else 1e-12
55 cy = c*(np.cos(ar) - np.cos(br)*np.cos(gr)) / denom
56 cz2 = c**2 - cx**2 - cy**2
57 cz = math.sqrt(max(cz2, 0.0))
58 c_vec = np.array([cx, cy, cz])
59 return np.column_stack((a_vec, b_vec, c_vec)) # 列: a,b,c
60
61def primitive_from_centering(A_conv, lattice_tag):
62 """
63 従来の格子基底と格子タイプから原始格子基底を計算します。
64
65 与えられた従来の格子基底に対し、指定されたセンタリングタイプ（P, I, Fなど）に
66 対応する変換行列を適用して原始格子基底を導出します。
67
68 :param A_conv: 従来の3x3格子基底行列（列が基底ベクトル）
69 :type A_conv: numpy.ndarray
70 :param lattice_tag: 格子タイプを示す文字列 (例: 'P', 'I', 'F', 'SC', 'BCC', 'FCC', 'A', 'B', 'C')
71 :type lattice_tag: str
72 :returns: 原始3x3格子基底行列（列が原始基底ベクトル）
73 :rtype: numpy.ndarray
74 :raises ValueError: 未知の格子タイプまたは未サポートの格子タイプが指定された場合
75 """
76 tag = lattice_tag.strip().upper()
77 alias = {'SC':'P','P':'P','BCC':'I','I':'I','FCC':'F','F':'F','A':'A','B':'B','C':'C'}
78 if tag not in alias: raise ValueError(f"Unknown lattice '{lattice_tag}'")
79 tag = alias[tag]
80 if tag == 'P':
81 P = np.eye(3)
82 elif tag == 'I':
83 P = np.array([[-0.5, 0.5, 0.5],
84 [ 0.5, -0.5, 0.5],
85 [ 0.5, 0.5, -0.5]])
86 elif tag == 'F':
87 P = np.array([[0.0, 0.5, 0.5],
88 [0.5, 0.0, 0.5],
89 [0.5, 0.5, 0.0]])
90 elif tag == 'R':
91 P = np.array([[ 2/3, -1/3, -1/3],
92 [ 1/3, 1/3, -2/3],
93 [ 1/3, 1/3, 1/3]])
94 elif tag == 'A':
95 P = np.array([[1.0, 0.0, 0.0],
96 [0.0, 0.5, 0.5],
97 [0.0, -0.5, 0.5]])
98 elif tag == 'B':
99 P = np.array([[0.5, 0.0, 0.5],
100 [0.0, 1.0, 0.0],
101 [0.5, 0.0, -0.5]])
102 elif tag == 'C':
103 P = np.array([[ 0.5, 0.5, 0.0],
104 [ 0.5, -0.5, 0.0],
105 [ 0.0, 0.0, 1.0]])
106 else:
107 raise ValueError(f"Unsupported lattice '{lattice_tag}'")
108 return A_conv @ P # 列が原始基底
109
110def order_polygon_vertices_on_plane(P):
111 """
112 平面上にある多角形の頂点を反時計回りに並べ替えるインデックスを返します。
113
114 頂点の重心を計算し、SVD（特異値分解）を用いて平面の法線ベクトルを決定します。
115 その後、基準ベクトルと法線に垂直なベクトルを生成し、これらの座標系で各頂点の角度を計算してソートします。
116 この順序はMatplotlibで多角形を描画する際に重要です。
117
118 :param P: 多角形の頂点群を表すNx3のNumPy配列。
119 :type P: numpy.ndarray
120 :returns: 頂点を反時計回りに並べ替えるためのインデックスの配列。
121 :rtype: numpy.ndarray
122 """
123 C = P.mean(axis=0); Q = P - C
124 if len(P) <= 2: return np.arange(len(P))
125 _, _, Vt = np.linalg.svd(Q, full_matrices=False)
126 n = Vt[-1]
127 ref = Q[0]
128 if np.linalg.norm(np.cross(ref, n)) < 1e-12:
129 ref = Q[1] if len(Q)>1 else np.array([1.0,0.0,0.0])
130 u = ref - np.dot(ref, n)*n
131 if np.linalg.norm(u) < 1e-12:
132 tmp = np.array([1.0,0.0,0.0]);
133 if abs(np.dot(tmp,n)) > 0.9: tmp = np.array([0.0,1.0,0.0])
134 u = tmp - np.dot(tmp, n)*n
135 u /= np.linalg.norm(u)
136 v = np.cross(n, u)
137 ang = np.arctan2(Q @ v, Q @ u)
138 return np.argsort(ang)
139
140def reciprocal_points_from_basis(B, n):
141 """
142 与えられた基底ベクトルとシェル数に基づいて、格子点または逆格子点の配列を生成します。
143
144 各基底ベクトルの整数倍（-nからnまで）の組み合わせで、中心の点を含む広範囲の点を生成します。
145 これはVoronoi図の計算に必要な点群を提供します。
146
147 :param B: 3x3の基底ベクトル行列（列が基底ベクトル）
148 :type B: numpy.ndarray
149 :param n: 各基底ベクトルに沿って考慮するシェル数（例: n=1なら-1,0,1の係数を使用）
150 :type n: int
151 :returns: 生成された点群のNx3のNumPy配列
152 :rtype: numpy.ndarray
153 """
154 pts = []
155 b1,b2,b3 = B[:,0], B[:,1], B[:,2]
156 for i in range(-n,n+1):
157 for j in range(-n,n+1):
158 for k in range(-n,n+1):
159 pts.append(i*b1 + j*b2 + k*b3)
160 return np.array(pts)
161
162def get_cell_faces_and_neighbors(vor, point_index):
163 """
164 Voronoi図の中心点に関連する面の頂点インデックスと隣接点のインデックスを取得します。
165
166 Voronoiオブジェクトのridge_pointsとridge_verticesを走査し、
167 指定された点に関連する有限な面（-1を含まない、3頂点以上）を抽出します。
168
169 :param vor: scipy.spatial.Voronoiオブジェクト
170 :type vor: scipy.spatial.Voronoi
171 :param point_index: 中心となる点の、Voronoi計算に用いた点群におけるインデックス
172 :type point_index: int
173 :returns: (faces, neighbors) のタプル。
174 facesは各面の頂点インデックスのリストのリスト、
175 neighborsは各面に対応する隣接点のインデックスのリスト。
176 :rtype: tuple[list[numpy.ndarray], list[int]]
177 """
178 faces, neighbors = [], []
179 for (p,q), v_idx in zip(vor.ridge_points, vor.ridge_vertices):
180 if point_index not in (p,q): continue
181 if -1 in v_idx or len(v_idx) < 3: continue
182 faces.append(np.array(v_idx))
183 neighbors.append(q if p==point_index else p)
184 return faces, neighbors
185
186def cell_polygons(vor, point_index):
187 """
188 Voronoi図の中心点に関連する多角形（セル面）と隣接点のリストを生成します。
189
190 `get_cell_faces_and_neighbors` を利用して面の頂点インデックスを取得し、
191 `order_polygon_vertices_on_plane` で頂点を反時計回りに順序付けて多角形オブジェクトを構築します。
192
193 :param vor: scipy.spatial.Voronoiオブジェクト
194 :type vor: scipy.spatial.Voronoi
195 :param point_index: 中心となる点の、Voronoi計算に用いた点群におけるインデックス
196 :type point_index: int
197 :returns: (polygons, neighbors) のタプル。
198 polygonsは各面の頂点座標のNx3 NumPy配列のリスト、
199 neighborsは各面に対応する隣接点のインデックスのリスト。
200 :rtype: tuple[list[numpy.ndarray], list[int]]
201 """
202 polys = []
203 faces, neighbors = get_cell_faces_and_neighbors(vor, point_index)
204 for f in faces:
205 poly = vor.vertices[f]
206 poly = poly[order_polygon_vertices_on_plane(poly)]
207 polys.append(poly)
208 return polys, neighbors
209
210# ================== SeeK-path（キャッシュ付き） ==================
211def _key_from_Aprim(A_prim, symprec):
212 """
213 原始格子基底行列と対称操作の許容誤差からキャッシュキーを生成します。
214
215 `functools.lru_cache` で使用するためのキーです。
216 浮動小数点数の比較を安定させるため、行列の要素は指定された精度で丸められます。
217
218 :param A_prim: 原始3x3格子基底行列（列が基底ベクトル）
219 :type A_prim: numpy.ndarray
220 :param symprec: 対称操作の許容誤差
221 :type symprec: float
222 :returns: キャッシュキーとして使用されるタプル
223 :rtype: tuple
224 """
225 # 行が a1,a2,a3 のフラット（数値は丸めてキャッシュ安定化）
226 Arows = A_prim.T
227 return (float(symprec),) + tuple(np.round(Arows.ravel(), 12))
228
229@functools.lru_cache(maxsize=128)
230def _seekpath_cached(key):
231 """
232 seekpath.get_path を呼び出し、結果をキャッシュする内部関数です。
233
234 `functools.lru_cache` デコレータにより、同じ入力（キャッシュキー）に対しては
235 実際に `seekpath` を呼び出さずに、キャッシュされた結果を返します。
236 これにより、インタラクティブな描画での再計算を高速化します。
237
238 :param key: `_key_from_Aprim` で生成されたキャッシュキー
239 :type key: tuple
240 :returns: (Bstd, kpts_abs, path) のタプル。
241 Bstdは標準逆格子基底行列、
242 kpts_absは高対称点の絶対座標辞書、
243 pathは高対称点パスのリスト。
244 :rtype: tuple[numpy.ndarray, dict[str, numpy.ndarray], list[tuple[str, str]]]
245 """
246 symprec = key[0]
247 Arows = np.array(key[1:]).reshape(3,3)
248 import seekpath
249 cell_tuple = (Arows, np.array([[0.0,0.0,0.0]]), [1])
250 res = seekpath.get_path(cell_tuple, with_time_reversal=True, symprec=symprec)
251 B_rows = np.array(res['reciprocal_primitive_lattice'])
252 Bstd = B_rows.T # 列が b1,b2,b3（2π含む）
253 kpts_abs = {lbl: Bstd @ np.array(frac) for lbl, frac in res['point_coords'].items()}
254 path = res['path']
255 return Bstd, kpts_abs, path
256
257def bz_and_kpath_from_seekpath_cached(A_prim, symprec=1e-7):
258 """
259 `seekpath` を用いて、第一ブリルアンゾーンの標準逆格子基底、高対称点、およびパスを取得します。
260
261 `_seekpath_cached` を通じて `functools.lru_cache` によるキャッシュを利用するため、
262 同じ原始格子基底に対しては `seekpath` の計算がスキップされ、高速に結果を返します。
263
264 :param A_prim: 原始3x3格子基底行列（列が基底ベクトル）
265 :type A_prim: numpy.ndarray
266 :param symprec: 対称操作の許容誤差
267 :type symprec: float
268 :returns: (Bstd, kpts_abs, path) のタプル。
269 Bstdは標準逆格子基底行列、
270 kpts_absは高対称点の絶対座標辞書、
271 pathは高対称点パスのリスト。
272 :rtype: tuple[numpy.ndarray, dict[str, numpy.ndarray], list[tuple[str, str]]]
273 """
274 key = _key_from_Aprim(A_prim, symprec)
275 return _seekpath_cached(key)
276
277def pretty_label(lbl):
278 """
279 高対称点のラベルを整形します（例: "GAMMA" を "Γ" に変換）。
280
281 :param lbl: 整形前の高対称点ラベル文字列
282 :type lbl: str
283 :returns: 整形された高対称点ラベル文字列
284 :rtype: str
285 """
286 return lbl.replace("GAMMA", "Γ")
287
288# ================== 着色＆透明度ユーティリティ ==================
289def color_checker(i,j,k):
290 """
291 三次元インデックス (i,j,k) に基づいてチェッカーボード状の色を返します。
292
293 各インデックスの偶奇性から0から7までのインデックスを生成し、
294 定義されたパレットから色を選択します。
295
296 :param i: 第1軸方向のインデックス
297 :type i: int
298 :param j: 第2軸方向のインデックス
299 :type j: int
300 :param k: 第3軸方向のインデックス
301 :type k: int
302 :returns: HEX形式の色コード文字列 (例: '#1f77b4')
303 :rtype: str
304 """
305 palette = ['#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd','#8c564b','#e377c2','#7f7f7f']
306 idx = ((i & 1) << 0) | ((j & 1) << 1) | ((k & 1) << 2) # 0..7
307 return palette[idx]
308
309def color_hash(i,j,k, cmap=plt.cm.tab20):
310 """
311 三次元インデックス (i,j,k) に基づいてハッシュ値から色を生成します。
312
313 与えられた整数インデックスから一意のハッシュ値を生成し、
314 そのハッシュ値を正規化して指定されたカラーマップから色を抽出します。
315
316 :param i: 第1軸方向のインデックス
317 :type i: int
318 :param j: 第2軸方向のインデックス
319 :type j: int
320 :param k: 第3軸方向のインデックス
321 :type k: int
322 :param cmap: 使用するMatplotlibカラーマップ
323 :type cmap: matplotlib.colors.Colormap
324 :returns: RGBA形式の色タプル
325 :rtype: tuple
326 """
327 h = (i*73856093 ^ j*19349663 ^ k*83492791) & 0xffffffff
328 return cmap((h % 256)/255.0)
329
330def color_distance(t, rmax, cmap=plt.cm.viridis):
331 """
332 原点からの距離に基づいてグラデーション色を返します。
333
334 点の位置ベクトル `t` の原点からの距離を計算し、`rmax` で正規化して
335 指定されたカラーマップ `cmap` から色を抽出します。
336
337 :param t: 点の3D座標ベクトル
338 :type t: numpy.ndarray
339 :param rmax: 色のグラデーションの最大半径。これを超えると色が飽和します。
340 Noneまたは0以下の場合は常に最大値の色を返します。
341 :type rmax: float or None
342 :param cmap: 使用するMatplotlibカラーマップ
343 :type cmap: matplotlib.colors.Colormap
344 :returns: RGBA形式の色タプル
345 :rtype: tuple
346 """
347 r = np.linalg.norm(t); u = 0.0 if rmax is None or rmax<=0 else min(1.0, r/rmax)
348 return cmap(u)
349
350def alpha_fade(t, rmax, a_min=0.06, a_max=0.22, power=1.0):
351 """
352 原点からの距離に基づいて透明度をフェードさせる値を返します。
353
354 点の位置ベクトル `t` の原点からの距離が `rmax` に近づくにつれて、
355 透明度が `a_max` から `a_min` へとフェードアウトするように計算されます。
356 `power` パラメータでフェードの非線形性を調整できます。
357
358 :param t: 点の3D座標ベクトル
359 :type t: numpy.ndarray
360 :param rmax: 透明度フェードの最大半径。これを超えると透明度は `a_min` に近づきます。
361 Noneまたは0以下の場合は常に `a_max` を返します。
362 :type rmax: float or None
363 :param a_min: 最小透明度 (0.0=完全透明, 1.0=完全不透明)
364 :type a_min: float
365 :param a_max: 最大透明度 (0.0=完全透明, 1.0=完全不透明)
366 :type a_max: float
367 :param power: フェードの指数（例: 1.0=線形, 2.0=二乗）
368 :type power: float
369 :returns: 計算された透明度の値
370 :rtype: float
371 """
372 if rmax is None or rmax<=0: return a_max
373 r = np.linalg.norm(t); u = min(1.0, r/rmax); w = (1.0 - u)**power
374 return a_min + (a_max - a_min)*w
375
376# ================== 描画マネージャ（差分更新） ==================
377class SceneManager:
378 """
379 Matplotlibの3Dシーンを管理し、描画オブジェクト（アーティスト）の追加、クリア、更新を行います。
380
381 特に、描画のたびに既存のオブジェクトを効率的に削除し、新しいオブジェクトを追加することで、
382 インタラクティブな描画でのちらつきを抑え、高速な更新を可能にします。
383 """
384 def __init__(self, args):
385 """
386 シーンマネージャを初期化します。
387
388 MatplotlibのFigureと3D Axesを作成し、以前の描画オブジェクトを保持するためのリストを準備します。
389
390 :param args: コマンドライン引数を格納するargparse.Namespaceオブジェクト
391 :type args: argparse.Namespace
392 """
393 self.args = args
394 self.fig = plt.figure(figsize=(10,9))
395 self.ax = self.fig.add_subplot(111, projection='3d')
396 self.artists = [] # 前回描画したアーティストを保持（差し替え用）
397
398 def clear_artists(self):
399 """
400 前回描画したすべてのMatplotlibアーティストをシーンから削除します。
401
402 これにより、Axes自体は維持しつつ、内容だけを更新することができます。
403 """
404 # Axesは維持し、前回のアーティストだけ remove()
405 for a in self.artists:
406 try: a.remove()
407 except Exception: pass
408 self.artists = []
409
410 def add_artist(self, *objs):
411 """
412 新しいMatplotlibアーティストを管理リストに追加します。
413
414 描画されたアーティストは、後で `clear_artists` メソッドで効率的に削除できるように
415 このリストに格納されます。リストやタプルの形式で複数のアーティストを渡すことも可能です。
416
417 :param objs: シーンに追加するMatplotlibアーティストオブジェクト（一つまたは複数）
418 :type objs: matplotlib.artist.Artist or list or tuple
419 """
420 for o in objs:
421 if isinstance(o, (list, tuple)):
422 self.artists.extend(o)
423 elif o is not None:
424 self.artists.append(o)
425
426 def draw_scene(self, space, A_prim, Bstd, n_shell,
427 tile_nx, tile_ny, tile_nz, tile_rmax,
428 show_outside_points, highlight_face_points,
429 draw_connectors_to_faces, draw_kpath,
430 color_mode='checker', fade_power=1.0,
431 kpts_abs=None, kpath=None,
432 label_fontsize=10):
433 """
434 第一ブリルアンゾーンまたはウィグナー・ザイツセル、およびそのタイリングされたビューを描画します。
435
436 空間タイプ（実空間または逆空間）に応じてVoronoi図を構築し、
437 中心セル、面生成点、周期タイル、高対称点とパス（逆空間の場合のみ）を描画します。
438 既存のアーティストをクリアし、新しいアーティストを追加する差分更新方式を採用しています。
439
440 :param space: 描画する空間タイプ ('reciprocal' または 'real')
441 :type space: str
442 :param A_prim: 原始3x3格子基底行列（実空間の基底）
443 :type A_prim: numpy.ndarray
444 :param Bstd: 標準逆格子3x3基底行列（逆空間の基底、2π因子を含む）
445 :type Bstd: numpy.ndarray or None
446 :param n_shell: Voronoi計算のために考慮する格子点/逆格子点のシェル数
447 :type n_shell: int
448 :param tile_nx: x軸方向のタイルの最大数（0で非表示）
449 :type tile_nx: int
450 :param tile_ny: y軸方向のタイルの最大数（0で非表示）
451 :type tile_ny: int
452 :param tile_nz: z軸方向のタイルの最大数（0で非表示）
453 :type tile_nz: int
454 :param tile_rmax: タイル描画の半径上限（中心から距離 r <= rmax のセルのみ描画）
455 Noneまたは0以下の場合、距離による制限なし
456 :type tile_rmax: float or None
457 :param show_outside_points: セル外部の格子点/逆格子点を表示するかどうか
458 :type show_outside_points: bool
459 :param highlight_face_points: セルの面を定義する点を強調表示するかどうか
460 :type highlight_face_points: bool
461 :param draw_connectors_to_faces: 面を定義する点と原点を結ぶ破線を描画するかどうか
462 :type draw_connectors_to_faces: bool
463 :param draw_kpath: 高対称点パスを描画するかどうか（reciprocal空間のみ有効）
464 :type draw_kpath: bool
465 :param color_mode: タイルの色付けモード ('checker', 'byindex', 'distance')
466 :type color_mode: str
467 :param fade_power: 距離に応じた透明度フェードの指数
468 :type fade_power: float
469 :param kpts_abs: 高対称点の絶対座標を格納した辞書 (reciprocal空間の場合のみ)
470 :type kpts_abs: dict[str, numpy.ndarray] or None
471 :param kpath: 高対称点パスのリスト (reciprocal空間の場合のみ)
472 :type kpath: list[tuple[str, str]] or None
473 :param label_fontsize: 高対称点ラベルのフォントサイズ
474 :type label_fontsize: int
475 """
476 self.clear_artists()
477 ax = self.ax
478
479 if space == 'reciprocal':
480 G = reciprocal_points_from_basis(Bstd, n_shell)
481 vor = Voronoi(G)
482 idx_center = np.argmin(np.linalg.norm(G, axis=1))
483 polys, neighbor_idx = cell_polygons(vor, idx_center)
484 tile_basis = Bstd
485 label_axes = (r'$k_x$', r'$k_y$', r'$k_z$')
486 outside_pts = G
487 face_color_main = 'lightskyblue'
488 origin_label = 'Γ (origin)'
489 else:
490 R = reciprocal_points_from_basis(A_prim, n_shell)
491 vor = Voronoi(R)
492 idx_center = np.argmin(np.linalg.norm(R, axis=1))
493 polys, neighbor_idx = cell_polygons(vor, idx_center)
494 tile_basis = A_prim
495 label_axes = ('x','y','z')
496 outside_pts = R
497 face_color_main = 'peachpuff'
498 origin_label = 'origin'
499
500 # 外部点群
501 if show_outside_points:
502 mask = np.ones(len(outside_pts), dtype=bool); mask[idx_center] = False
503 sc = ax.scatter(outside_pts[mask,0], outside_pts[mask,1], outside_pts[mask,2],
504 s=10, color='gray', alpha=0.45,
505 label=('Reciprocal points' if space=='reciprocal' else 'Lattice points'))
506 self.add_artist(sc)
507
508 # 中心セル（外周のみ）
509 for poly in polys:
510 col = Poly3DCollection([poly], alpha=0.32, facecolor=face_color_main, edgecolor='none')
511 ax.add_collection3d(col)
512 self.add_artist(col)
513 ring = np.vstack([poly, poly[0]])
514 ln, = ax.plot(ring[:,0], ring[:,1], ring[:,2], color='k', linewidth=1.8)
515 self.add_artist(ln)
516
517 # 面生成点と点線
518 centers = outside_pts[np.array(neighbor_idx)] if len(neighbor_idx)>0 else np.empty((0,3))
519 if highlight_face_points and len(centers)>0:
520 sc2 = ax.scatter(centers[:,0], centers[:,1], centers[:,2], s=28, color='tab:blue', alpha=0.9, label='Face-defining points')
521 self.add_artist(sc2)
522 if draw_connectors_to_faces:
523 lines = []
524 for p in centers:
525 ln, = ax.plot([0,p[0]],[0,p[1]],[0,p[2]], linestyle=':', linewidth=1.2, color='k', alpha=0.9)
526 lines.append(ln)
527 self.add_artist(lines)
528
529 # 周期タイル
530 if any(n>0 for n in (tile_nx, tile_ny, tile_nz)):
531 b1,b2,b3 = tile_basis[:,0], tile_basis[:,1], tile_basis[:,2]
532 for i in range(-tile_nx, tile_nx+1):
533 for j in range(-tile_ny, tile_ny+1):
534 for k in range(-tile_nz, tile_nz+1):
535 if i==0 and j==0 and k==0: continue
536 t = i*b1 + j*b2 + k*b3
537 if tile_rmax is not None and np.linalg.norm(t) > tile_rmax + 1e-12:
538 continue
539 # 色
540 if color_mode == 'checker':
541 fc = color_checker(i,j,k)
542 elif color_mode == 'distance':
543 fc = color_distance(t, tile_rmax)
544 else:
545 fc = color_hash(i,j,k)
546 # 透明度
547 a = alpha_fade(t, tile_rmax, a_min=0.06, a_max=0.22, power=fade_power)
548 # セル描画
549 cols, lns = [], []
550 for poly in polys:
551 P = poly + t
552 c = Poly3DCollection([P], alpha=a, facecolor=fc, edgecolor='none')
553 ax.add_collection3d(c); cols.append(c)
554 ring = np.vstack([P, P[0]])
555 ln, = ax.plot(ring[:,0], ring[:,1], ring[:,2], color='#666666', linewidth=0.6); lns.append(ln)
556 self.add_artist(cols); self.add_artist(lns)
557
558 # 高対称点と k-path（逆格子のみ）
559 if space=='reciprocal' and kpts_abs:
560 texts, pts, segs = [], [], []
561 for lbl, k in kpts_abs.items():
562 sct = ax.scatter(k[0],k[1],k[2], s=40, color='crimson'); pts.append(sct)
563 txt = ax.text(k[0],k[1],k[2], " "+pretty_label(lbl), fontsize=label_fontsize, color='crimson'); texts.append(txt)
564 self.add_artist(pts); self.add_artist(texts)
565 if draw_kpath and kpath:
566 for a_lbl,b_lbl in kpath:
567 ka,kb = kpts_abs[a_lbl], kpts_abs[b_lbl]
568 seg = np.vstack([ka,kb])
569 ln, = ax.plot(seg[:,0], seg[:,1], seg[:,2], color='crimson', linewidth=1.6, alpha=0.85)
570 segs.append(ln)
571 self.add_artist(segs)
572
573 # 原点
574 origin = ax.scatter([0],[0],[0], color='k', s=50, label=origin_label)
575 self.add_artist(origin)
576
577 # 軸や凡例は1回だけ設定（差分更新では再設定しない）
578 basis = (Bstd if space=='reciprocal' else A_prim)
579 scale = 1.35 * np.linalg.norm(basis[:,0])
580 ax.set_xlim(-scale, scale); ax.set_ylim(-scale, scale); ax.set_zlim(-scale, scale)
581
582 if self.args.plot_limit is not None:
583 L = self.args.plot_limit
584 ax.set_xlim(-L, L)
585 ax.set_ylim(-L, L)
586 ax.set_zlim(-L, L)
587 ax.set_box_aspect([1,1,1])
588 ax.set_xlabel(label_axes[0]); ax.set_ylabel(label_axes[1]); ax.set_zlabel(label_axes[2])
589 ax.set_title('First Brillouin Zone (tiled)' if space=='reciprocal' else 'Wigner–Seitz cell (tiled)')
590 # 凡例は毎回作り直す（重複回避のため一旦消してから）
591 leg = ax.legend(loc='upper right'); self.add_artist(leg)
592
593 self.fig.tight_layout()
594 self.fig.canvas.draw_idle()
595
596# ================== 実行系 ==================
597def plot_once(args):
598 """
599 指定されたパラメータに基づいて、第一ブリルアンゾーンまたはウィグナー・ザイツセルを一度だけ描画します。
600
601 結晶学的パラメータから格子基底を構築し、`SceneManager` を使用してシーンを描画し、
602 結果をMatplotlibウィンドウに表示します。
603
604 :param args: コマンドライン引数を格納するargparse.Namespaceオブジェクト
605 :type args: argparse.Namespace
606 :raises SystemExit: `seekpath` の実行に失敗した場合
607 """
608 A_conv = build_conventional_basis(args.a, args.b, args.c, args.alpha, args.beta, args.gamma)
609 A_prim = primitive_from_centering(A_conv, args.lattice)
610
611 mgr = SceneManager(args)
612 if args.space == 'reciprocal':
613 try:
614 Bstd, kpts_abs, kpath = bz_and_kpath_from_seekpath_cached(A_prim)
615 except RuntimeError as e:
616 print(e, file=sys.stderr); sys.exit(1)
617 mgr.draw_scene('reciprocal', A_prim, Bstd, args.n_shell,
618 args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
619 (not args.no_outside_points), (not args.no_highlight_face_points),
620 (not args.no_connectors), (not args.no_kpath),
621 color_mode=args.tile_color_mode, fade_power=args.fade_power,
622 kpts_abs=kpts_abs, kpath=kpath, label_fontsize=args.label_fontsize)
623 else:
624 mgr.draw_scene('real', A_prim, None, args.n_shell,
625 args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
626 (not args.no_outside_points), (not args.no_highlight_face_points),
627 (not args.no_connectors), False,
628 color_mode=args.tile_color_mode, fade_power=args.fade_power,
629 kpts_abs=None, kpath=None, label_fontsize=args.label_fontsize)
630 plt.show()
631
632def plot_interactive_apply(args):
633 """
634 スライダーと「Apply」ボタンを備えたインタラクティブな描画インターフェースを提供します。
635
636 `a, b, c` の格子定数をスライダーで調整でき、「Apply」ボタンをクリックしたときにのみ
637 シーンが再描画されます。`seekpath` の結果はキャッシュされるため、
638 再描画時の計算コストが削減され、高速な更新が可能です。
639
640 :param args: コマンドライン引数を格納するargparse.Namespaceオブジェクト
641 :type args: argparse.Namespace
642 """
643 # 初期値
644 A_conv = build_conventional_basis(args.a, args.b, args.c, args.alpha, args.beta, args.gamma)
645 A_prim = primitive_from_centering(A_conv, args.lattice)
646 mgr = SceneManager(args)
647
648 def do_draw(a, b, c):
649 A_conv = build_conventional_basis(a, b, c, args.alpha, args.beta, args.gamma)
650 A_prim = primitive_from_centering(A_conv, args.lattice)
651 if args.space == 'reciprocal':
652 Bstd, kpts_abs, kpath = bz_and_kpath_from_seekpath_cached(A_prim) # キャッシュで高速化
653 mgr.draw_scene('reciprocal', A_prim, Bstd, args.n_shell,
654 args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
655 (not args.no_outside_points), (not args.no_highlight_face_points),
656 (not args.no_connectors), (not args.no_kpath),
657 color_mode=args.tile_color_mode, fade_power=args.fade_power,
658 kpts_abs=kpts_abs, kpath=kpath, label_fontsize=args.label_fontsize)
659 else:
660 mgr.draw_scene('real', A_prim, None, args.n_shell,
661 args.tile_nx, args.tile_ny, args.tile_nz, args.tile_rmax,
662 (not args.no_outside_points), (not args.no_highlight_face_points),
663 (not args.no_connectors), False,
664 color_mode=args.tile_color_mode, fade_power=args.fade_power,
665 kpts_abs=None, kpath=None, label_fontsize=args.label_fontsize)
666
667 # 初回描画
668 do_draw(args.a, args.b, args.c)
669
670 # UI: スライダ + Apply ボタン（押したときのみ再描画）
671 plt.subplots_adjust(bottom=0.20)
672 ax_a = plt.axes([0.15, 0.08, 0.65, 0.03])
673 ax_b = plt.axes([0.15, 0.05, 0.65, 0.03])
674 ax_c = plt.axes([0.15, 0.02, 0.65, 0.03])
675 s_a = Slider(ax_a, 'a', 0.5*args.a, 2.0*args.a, valinit=args.a)
676 s_b = Slider(ax_b, 'b', 0.5*args.b, 2.0*args.b, valinit=args.b)
677 s_c = Slider(ax_c, 'c', 0.5*args.c, 2.0*args.c, valinit=args.c)
678
679 ax_btn = plt.axes([0.83, 0.02, 0.12, 0.09])
680 btn = Button(ax_btn, 'Apply')
681
682 busy = {'flag': False}
683 def on_apply(event):
684 if busy['flag']: # 処理中は無視
685 return
686 busy['flag'] = True
687 # ボタン表示を Redrawing... に変更
688 btn.label.set_text('Redrawing...')
689 btn.ax.set_facecolor('#dddddd')
690 mgr.fig.canvas.draw_idle()
691
692 try:
693 do_draw(s_a.val, s_b.val, s_c.val)
694 finally:
695 # ボタン表示を戻す
696 btn.label.set_text('Apply')
697 btn.ax.set_facecolor('#f0f0f0')
698 mgr.fig.canvas.draw_idle()
699 busy['flag'] = False
700
701 btn.on_clicked(on_apply)
702 plt.show()
703
704# ================== CLI ==================
705def main():
706 """
707 プログラムのエントリーポイントです。
708
709 コマンドライン引数を解析し、ユーザーの選択に応じて
710 第一ブリルアンゾーンまたはウィグナー・ザイツセルの描画機能
711 （一度きりの描画またはインタラクティブな描画）を呼び出します。
712 """
713 ap = argparse.ArgumentParser(description="Voronoi 多面体（第一BZ / Wigner–Seitz）: キャッシュ・差分更新・Applyボタン対応")
714 ap.add_argument('--space', choices=['reciprocal','real'], default='reciprocal',
715 help="reciprocal: 第一BZ（SeeK-pathで高対称点） / real: 実空間 Wigner–Seitz")
716 ap.add_argument('--lattice', default='P', help="格子タイプ: P|SC, I|BCC, F|FCC, A, B, C")
717 ap.add_argument('--a', type=float, default=1.0); ap.add_argument('--b', type=float, default=1.0); ap.add_argument('--c', type=float, default=1.0)
718 ap.add_argument('--alpha', type=float, default=90.0); ap.add_argument('--beta', type=float, default=90.0); ap.add_argument('--gamma', type=float, default=90.0)
719 ap.add_argument('--n-shell', type=int, default=3, help='Voronoi 計算用に並べる格子点の範囲')
720 ap.add_argument('--tile-nx', type=int, default=0); ap.add_argument('--tile-ny', type=int, default=0); ap.add_argument('--tile-nz', type=int, default=0)
721 ap.add_argument('--tile-rmax', type=float, default=None,
722 help='タイル描画の半径上限（中心から距離 r ≤ rmax のセルのみ）')
723 ap.add_argument('--plot-limit', type=float, default=None,
724 help='軸の表示範囲を [-L, L] に強制的に制限する（内部座標、L=1.0で[-2pi/a, 2pi/a]に対応）')
725 ap.add_argument('--tile-color-mode', choices=['checker','byindex','distance'], default='checker',
726 help="タイル色: checker=チェッカー, byindex=ハッシュ, distance=距離グラデーション")
727 ap.add_argument('--fade-power', type=float, default=1.0, help='距離フェードの指数（1=線形, 2=二乗など）')
728 ap.add_argument('--no-outside-points', action='store_true')
729 ap.add_argument('--no-highlight-face-points', action='store_true')
730 ap.add_argument('--no-connectors', action='store_true')
731 ap.add_argument('--no-kpath', action='store_true')
732 ap.add_argument('--label-fontsize', type=int, default=10)
733 ap.add_argument('--interactive-apply', action='store_true', help='スライダ + Applyボタンで再描画')
734 ap.add_argument('--save', default=None, help='（一括描画時のみ）画像保存パス')
735
736 args = ap.parse_args()
737 if args.interactive_apply:
738 plot_interactive_apply(args)
739 else:
740 plot_once(args)
741
742if __name__ == "__main__":
743 main()