draw_bravais_lattice.py ダウンロード/コピー
draw_bravais_lattice.py をダウンロード
draw_bravais_lattice.py
draw_bravais_lattice.py
1"""
2ブラベー格子の単位格子と格子点を3Dで可視化するスクリプト。
3
4このスクリプトは、指定された格子定数と結晶系に基づいて、
5ブラベー格子の単位格子、格子点、格子ベクトル、およびオプションで補助線や基本格子を3Dで描画します。
6Matplotlibを使用して3Dプロットを生成し、結果を画像ファイルとして保存し、表示します。
7
8関連リンク: :doc:`draw_bravais_lattice_usage`
9"""
10import sys
11import numpy as np
12import matplotlib.pyplot as plt
13from mpl_toolkits.mplot3d import Axes3D
14from matplotlib.patches import FancyArrowPatch
15from mpl_toolkits.mplot3d import Axes3D, proj3d
16
17# 始点の初期値
18elev = 5.90
19azim = -67.55
20
21draw_lattice_vectors = True
22draw_support_lines = True
23draw_primitive_cell = False
24
25# ブラべー格子
26lattice_vector_color = 'blue'
27axis_arrow_color = 'blue'
28axis_label_font_size = 24
29# 格子点
30point_color = 'red'
31point_size = 500
32
33#cell = 'SC'
34#cell = 'ST'
35#cell = 'SO'
36#cell = 'SR'
37#cell = 'SH'
38#cell = 'FC'
39#cell = 'FO'
40cell = 'BC'
41#cell = 'BT'
42#cell = 'BO'
43#cell = 'CO'
44#cell = 'SM'
45#cell = 'STri'
46#cell = 'CM'
47
48
49class Arrow3D(FancyArrowPatch):
50 """
51 MatplotlibのFancyArrowPatchを継承し、3D空間に矢印を描画するためのクラス。
52
53 このクラスは、3D座標の矢印を2Dのスクリーン座標に投影し、
54 Matplotlibの3D軸上に適切に描画することを可能にします。
55 """
56 def __init__(self, xs, ys, zs, *args, **kwargs):
57 """
58 Arrow3Dオブジェクトを初期化します。
59
60 :param xs: 矢印のX座標 (開始点, 終了点)。
61 :type xs: list[float]
62 :param ys: 矢印のY座標 (開始点, 終了点)。
63 :type ys: list[float]
64 :param zs: 矢印のZ座標 (開始点, 終了点)。
65 :type zs: list[float]
66 :param args: FancyArrowPatchに渡される追加の位置引数。
67 :param kwargs: FancyArrowPatchに渡される追加のキーワード引数。
68 """
69 super().__init__((0, 0), (0, 0), *args, **kwargs)
70 self._verts3d = xs, ys, zs
71
72 def draw(self, renderer):
73 """
74 矢印をレンダラーに描画します。
75
76 3D座標を2Dスクリーン座標に変換し、変換された座標に基づいて矢印を描画します。
77
78 :param renderer: 描画に使用されるレンダラーオブジェクト。
79 :type renderer: matplotlib.backend_bases.RendererBase
80 """
81 xs3d, ys3d, zs3d = self._verts3d
82 proj = self.axes.get_proj()
83 xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
84 self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
85 super().draw(renderer)
86
87 def do_3d_projection(self, renderer=None):
88 """
89 3D投影を計算し、描画順序のためのZ深度を返します。
90
91 矢印の3D座標を2Dスクリーン座標に変換し、最小のZ深度を返して、
92 Matplotlibの3Dレンダリングにおける描画順序を適切に処理します。
93
94 :param renderer: 描画に使用されるレンダラーオブジェクト（オプション）。
95 :type renderer: matplotlib.backend_bases.RendererBase or None
96 :returns: 変換されたZ座標の最小値。
97 :rtype: float
98 """
99 xs3d, ys3d, zs3d = self._verts3d
100 proj = self.axes.get_proj()
101 xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, proj)
102 self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
103 return np.min(zs)
104
105
106def draw_vector(ax, vec, color, label, fontsize):
107 """
108 指定された3Dベクトルを矢印としてMatplotlibの3D軸に描画します。
109
110 矢印は原点(0,0,0)から指定されたベクトルまで描画され、
111 ベクトルの先端にはラベルが表示されます。
112
113 :param ax: 矢印を描画する3D軸オブジェクト。
114 :type ax: mpl_toolkits.mplot3d.axes3d.Axes3D
115 :param vec: 矢印の終点となる3Dベクトル。
116 :type vec: numpy.ndarray
117 :param color: 矢印とラベルの色。
118 :type color: str
119 :param label: 矢印のラベル文字列。
120 :type label: str
121 :param fontsize: ラベルのフォントサイズ。
122 :type fontsize: int
123 :returns: なし
124 """
125 arrow = Arrow3D([0, vec[0]], [0, vec[1]], [0, vec[2]],
126 mutation_scale=20, lw=2, arrowstyle="-|>", color=color)
127 ax.add_artist(arrow)
128 ax.text(vec[0] * 1.05, vec[1] * 1.05, vec[2] * 1.05, label, fontsize=fontsize, color=color)
129
130
131def lattice_vectors(a, b, c, alpha, beta, gamma):
132 """
133 単位格子の格子定数から3つの基本格子ベクトルを計算します。
134
135 計算は、a軸をx軸に沿わせ、b軸をxy平面に置き、c軸をz軸正方向に向ける形で
136 座標系を設定して行われます。
137
138 :param a: 格子定数aの長さ。
139 :type a: float
140 :param b: 格子定数bの長さ。
141 :type b: float
142 :param c: 格子定数cの長さ。
143 :type c: float
144 :param alpha: 角度α（bとcの間の角度、度数）。
145 :type alpha: float
146 :param beta: 角度β（aとcの間の角度、度数）。
147 :type beta: float
148 :param gamma: 角度γ（aとbの間の角度、度数）。
149 :type gamma: float
150 :returns: 3つの基本格子ベクトル (va, vb, vc)。
151 :rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]
152 """
153 alpha_r, beta_r, gamma_r = np.radians([alpha, beta, gamma])
154 va = np.array([a, 0, 0])
155 vb = np.array([b * np.cos(gamma_r), b * np.sin(gamma_r), 0])
156 cx = c * np.cos(beta_r)
157 cy = c * (np.cos(alpha_r) - np.cos(beta_r) * np.cos(gamma_r)) / np.sin(gamma_r)
158 cz = np.sqrt(c**2 - cx**2 - cy**2)
159 vc = np.array([cx, cy, cz])
160 return va, vb, vc
161
162
163def set_equal_aspect(ax, points):
164 """
165 Matplotlibの3Dプロットの軸のアスペクト比を等しく設定します。
166
167 描画される全ての点の範囲に基づいて、x, y, z軸の表示範囲を調整し、
168 立体が歪まないように真のアスペクト比を維持します。
169
170 :param ax: 軸のアスペクト比を設定する3D軸オブジェクト。
171 :type ax: mpl_toolkits.mplot3d.axes3d.Axes3D
172 :param points: 描画される全ての点のNumpy配列。
173 :type points: numpy.ndarray
174 :returns: なし
175 """
176 xlim = [np.min(points[:, 0]), np.max(points[:, 0])]
177 ylim = [np.min(points[:, 1]), np.max(points[:, 1])]
178 zlim = [np.min(points[:, 2]), np.max(points[:, 2])]
179 max_range = max(
180 xlim[1] - xlim[0],
181 ylim[1] - ylim[0],
182 zlim[1] - zlim[0]
183 ) / 2
184
185 mid_x = np.mean(xlim)
186 mid_y = np.mean(ylim)
187 mid_z = np.mean(zlim)
188
189 ax.set_xlim(mid_x - max_range, mid_x + max_range)
190 ax.set_ylim(mid_y - max_range, mid_y + max_range)
191 ax.set_zlim(mid_z - max_range, mid_z + max_range)
192 ax.set_box_aspect([1, 1, 1])
193
194
195def draw_unit_cell_with_lattice(
196 a, b, c, alpha, beta, gamma,
197 lattice_points=None,
198 dashed_lines=None,
199 basis_vectors=None
200):
201 """
202 指定された格子定数と格子点を用いて、ブラベー単位格子を3Dで描画します。
203
204 この関数は、単位格子の辺、頂点、格子ベクトル、格子点、オプションで補助線、
205 および基本格子を描画し、結果を画像ファイルとして保存した後、画面に表示します。
206 ユーザーがプロットを操作すると、現在の視点角度 (elev, azim) がコンソールに出力されます。
207
208 :param a: 格子定数aの長さ。
209 :type a: float
210 :param b: 格子定数bの長さ。
211 :type b: float
212 :param c: 格子定数cの長さ。
213 :type c: float
214 :param alpha: 角度α（bとcの間の角度、度数）。
215 :type alpha: float
216 :param beta: 角度β（aとcの間の角度、度数）。
217 :type beta: float
218 :param gamma: 角度γ（aとbの間の角度、度数）。
219 :type gamma: float
220 :param lattice_points: 格子点を表す係数のリスト。各要素は `[u, v, w]` 形式で `u*va + v*vb + w*vc` と計算されます。
221 Noneの場合、格子点は描画されません。
222 :type lattice_points: list[list[float]] or None
223 :param dashed_lines: 補助線（破線）の始点と終点を表す係数のリスト。各要素は `([u1, v1, w1], [u2, v2, w2])` 形式です。
224 Noneの場合、補助線は描画されません。
225 :type dashed_lines: list[tuple[tuple[float, float, float], tuple[float, float, float]]] or None
226 :param basis_vectors: 基本格子ベクトルを表す係数のリスト。各要素は `[u, v, w]` 形式です。
227 Noneの場合、基本格子は描画されません。
228 :type basis_vectors: list[list[float]] or None
229 :returns: なし
230 """
231 va, vb, vc = lattice_vectors(a, b, c, alpha, beta, gamma)
232 origin = np.array([0, 0, 0])
233 corners = [
234 origin,
235 va,
236 vb,
237 vc,
238 va + vb,
239 va + vc,
240 vb + vc,
241 va + vb + vc
242 ]
243 corners = np.array(corners)
244
245 edges = [
246 (0, 1), (0, 2), (0, 3),
247 (1, 4), (1, 5),
248 (2, 4), (2, 6),
249 (3, 5), (3, 6),
250 (4, 7), (5, 7), (6, 7)
251 ]
252
253 fig = plt.figure(figsize=(8, 6))
254 ax = fig.add_subplot(111, projection='3d')
255 ax.view_init(elev=elev, azim=azim)
256
257 def on_draw(event):
258 elev_now = ax.elev
259 azim_now = ax.azim
260 print(f"elev: {elev_now:.2f}, azim: {azim_now:.2f}")
261
262 fig.canvas.mpl_connect('draw_event', on_draw)
263
264 # 単位格子の辺
265 for i, j in edges:
266 ax.plot(*zip(corners[i], corners[j]), color='black', linewidth=0.5)
267
268 # 単位格子の頂点
269 ax.scatter(corners[:, 0], corners[:, 1], corners[:, 2], color='black', s=20)
270
271 # 格子ベクトル
272 if draw_lattice_vectors:
273 for vec, label in zip([va, vb, vc], ['a', 'b', 'c']):
274 draw_vector(ax, vec, lattice_vector_color, label, axis_label_font_size)
275
276 # 格子点の描画
277 real_points = []
278 if lattice_points:
279 real_points = [p[0] * va + p[1] * vb + p[2] * vc for p in lattice_points]
280 real_points = np.array(real_points)
281 ax.scatter(real_points[:, 0], real_points[:, 1], real_points[:, 2],
282 color=point_color, s=point_size)
283
284 # 補助線（鎖線）
285 if draw_support_lines and dashed_lines:
286 for p1, p2 in dashed_lines:
287 pt1 = p1[0] * va + p1[1] * vb + p1[2] * vc
288 pt2 = p2[0] * va + p2[1] * vb + p2[2] * vc
289 ax.plot(*zip(pt1, pt2), color='gray', linestyle='dashed', linewidth=0.8)
290
291 # 基本格子
292 if draw_primitive_cell and basis_vectors:
293 for vec in basis_vectors:
294 v = vec[0] * va + vec[1] * vb + vec[2] * vc
295 ax.quiver(0, 0, 0, *v, color='green', linewidth=1.5, arrow_length_ratio=0.2)
296
297 if len(basis_vectors) == 3:
298 bv1 = basis_vectors[0][0] * va + basis_vectors[0][1] * vb + basis_vectors[0][2] * vc
299 bv2 = basis_vectors[1][0] * va + basis_vectors[1][1] * vb + basis_vectors[1][2] * vc
300 bv3 = basis_vectors[2][0] * va + basis_vectors[2][1] * vb + basis_vectors[2][2] * vc
301
302 base_points = [
303 origin,
304 bv1,
305 bv2,
306 bv3,
307 bv1 + bv2,
308 bv1 + bv3,
309 bv2 + bv3,
310 bv1 + bv2 + bv3
311 ]
312
313 base_edges = [
314 (0, 1), (0, 2), (0, 3),
315 (1, 4), (1, 5),
316 (2, 4), (2, 6),
317 (3, 5), (3, 6),
318 (4, 7), (5, 7), (6, 7)
319 ]
320
321 for i, j in base_edges:
322 ax.plot(*zip(base_points[i], base_points[j]),
323 color='lightgreen', linewidth=1.2)
324
325 for vec, label in zip([bv1, bv2, bv3], ["a'", "b'", "c'"]):
326 ax.text(*vec, label, fontsize=12, color='green')
327
328 # 軸や背景の調整
329 ax.set_xticks([]); ax.set_yticks([]); ax.set_zticks([])
330 ax.set_xticklabels([]); ax.set_yticklabels([]); ax.set_zticklabels([])
331 ax.set_xlabel(''); ax.set_ylabel(''); ax.set_zlabel('')
332 ax.grid(False)
333 ax.xaxis.pane.set_visible(False)
334 ax.yaxis.pane.set_visible(False)
335 ax.zaxis.pane.set_visible(False)
336 ax.xaxis.line.set_color((0, 0, 0, 0))
337 ax.yaxis.line.set_color((0, 0, 0, 0))
338 ax.zaxis.line.set_color((0, 0, 0, 0))
339 ax.set_facecolor((1, 1, 1, 0))
340
341 all_points = np.vstack([corners] + ([real_points] if lattice_points else []))
342 set_equal_aspect(ax, all_points)
343
344 plt.tight_layout()
345 plt.savefig("bravais_cell_with_basis_vectors.png",
346 dpi=300, bbox_inches='tight', transparent=True)
347 plt.show()
348
349
350def main():
351 """
352 プログラムのエントリポイント。コマンドライン引数を解析し、ブラベー単位格子の描画を実行します。
353
354 コマンドライン引数に基づいて、描画する結晶系と表示オプション（格子ベクトル、補助線、基本格子）を
355 設定し、`draw_unit_cell_with_lattice`関数を呼び出して単位格子を描画します。
356 引数が指定されない場合は、デフォルト値が使用されます。
357
358 使用方法:
359 python draw_bravais_lattice.py [cell_type] [draw_lattice_vectors] [draw_support_lines] [draw_primitive_cell]
360
361 - cell_type (str): 描画する結晶系 ('SC', 'ST', 'SO', 'SR', 'SH', 'FC', 'FO', 'BC', 'BT', 'BO', 'CO', 'SM', 'STri', 'CM'など)。
362 - draw_lattice_vectors (int): 格子ベクトルを描画するか (1: はい, 0: いいえ)。
363 - draw_support_lines (int): 補助線を描画するか (1: はい, 0: いいえ)。
364 - draw_primitive_cell (int): 基本格子を描画するか (1: はい, 0: いいえ)。
365
366 :returns: なし
367 """
368 global cell, draw_lattice_vectors, draw_support_lines, draw_primitive_cell
369
370 argv = sys.argv
371 nargs = len(argv)
372 if nargs > 1:
373 cell = argv[1]
374 if nargs > 2:
375 draw_lattice_vectors = int(argv[2])
376 if nargs > 3:
377 draw_support_lines = int(argv[3])
378 if nargs > 4:
379 draw_primitive_cell = int(argv[4])
380
381 print()
382 print(f"{cell=}")
383 print(f"{draw_lattice_vectors=}")
384 print(f"{draw_support_lines=}")
385 print(f"{draw_primitive_cell=}")
386
387 if cell[0] == 'S': # Simple (単純)
388 if cell[1] == 'C': # Simple Cubic (単純立方)
389 a, b, c, alpha, beta, gamma = 5.0, 5.0, 5.0, 90.0, 90.0, 90.0
390 elif cell[1:] == 'T': # Simple Tetragonal (単純正方)
391 a, b, c, alpha, beta, gamma = 5.0, 5.0, 7.0, 90.0, 90.0, 90.0
392 elif cell[1] == 'O': # Simple Orthorhombic (単純斜方)
393 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 90.0, 90.0, 90.0
394 elif cell[1] == 'H': # Simple Hexagonal (単純六方)
395 a, b, c, alpha, beta, gamma = 7.0, 7.0, 5.0, 90.0, 90.0, 120.0
396 elif cell[1] == 'R': # Simple Rhombohedral (単純菱面体)
397 a, b, c, alpha, beta, gamma = 5.0, 5.0, 5.0, 70.0, 70.0, 70.0
398 elif cell[1] == 'M': # Simple Monoclinic (単純単斜)
399 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 90.0, 110.0, 90.0
400 elif cell[1:] == 'Tri': # Simple Triclinic (単純三斜)
401 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 80.0, 110.0, 60.0
402
403 lattice_points = [
404 [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1],
405 [0, 1, 1], [1, 0, 1], [1, 1, 0], [1, 1, 1],
406 ]
407
408 dashed_lines = []
409 basis_vectors = []
410
411 elif cell[0] == 'F': # Face-centered (面心)
412 if cell[1] == 'C': # Face-centered Cubic (面心立方)
413 a, b, c, alpha, beta, gamma = 5.0, 5.0, 5.0, 90.0, 90.0, 90.0
414 elif cell[1] == 'O': # Face-centered Orthorhombic (面心斜方)
415 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 90.0, 90.0, 90.0
416
417 lattice_points = [
418 [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1],
419 [0, 1, 1], [1, 0, 1], [1, 1, 0], [1, 1, 1],
420 [0.0, 0.5, 0.5],
421 [0.5, 0.0, 0.5],
422 [0.5, 0.5, 0.0],
423 [1.0, 0.5, 0.5],
424 [0.5, 1.0, 0.5],
425 [0.5, 0.5, 1.0],
426 ]
427
428 dashed_lines = [
429 [(0, 0, 0), (1, 1, 0)],
430 [(0, 1, 0), (1, 0, 0)],
431 [(0, 0, 1), (1, 1, 1)],
432 [(0, 1, 1), (1, 0, 1)],
433 [(0, 0, 0), (1, 0, 1)],
434 [(0, 0, 1), (1, 0, 0)],
435 [(0, 1, 0), (1, 1, 1)],
436 [(0, 1, 1), (1, 1, 0)],
437 [(0, 0, 0), (0, 1, 1)],
438 [(0, 1, 0), (0, 0, 1)],
439 [(1, 0, 0), (1, 1, 1)],
440 [(1, 1, 0), (1, 0, 1)],
441 ]
442
443 basis_vectors = [
444 [-0.5, 0.5, 0.5],
445 [ 0.5, -0.5, 0.5],
446 [ 0.5, 0.5, -0.5]
447 ]
448
449 elif cell[0] == 'B': # Body-centered (体心)
450 if cell[1] == 'C': # Body-centered Cubic (体心立方)
451 a, b, c, alpha, beta, gamma = 5.0, 5.0, 5.0, 90.0, 90.0, 90.0
452 elif cell[1] == 'T': # Body-centered Tetragonal (体心正方)
453 a, b, c, alpha, beta, gamma = 5.0, 5.0, 7.0, 90.0, 90.0, 90.0
454 elif cell[1] == 'O': # Body-centered Orthorhombic (体心斜方)
455 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 90.0, 90.0, 90.0
456
457 lattice_points = [
458 [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1],
459 [0, 1, 1], [1, 0, 1], [1, 1, 0], [1, 1, 1],
460 [0.5, 0.5, 0.5]
461 ]
462
463 dashed_lines = [
464 [(0, 0, 0), (1, 1, 1)],
465 [(1, 0, 0), (0, 1, 1)],
466 [(0, 1, 0), (1, 0, 1)],
467 [(0, 0, 1), (1, 1, 0)],
468 [(0.5, 0.5, 0.5), (0, 0, 0)],
469 [(0.5, 0.5, 0.5), (1, 1, 1)]
470 ]
471
472 basis_vectors = [
473 [-0.5, 0.5, 0.5],
474 [ 0.5, -0.5, 0.5],
475 [ 0.5, 0.5, -0.5]
476 ]
477
478 elif cell[0] == 'C': # Base-centered (底心)
479 if cell[1] == 'O': # Base-centered Orthorhombic (底心斜方)
480 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 90.0, 90.0, 90.0
481 lattice_points = [
482 [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1],
483 [0, 1, 1], [1, 0, 1], [1, 1, 0], [1, 1, 1],
484 [0.5, 0.5, 0.0],
485 [0.5, 0.5, 1.0]
486 ]
487 dashed_lines = [
488 [(0, 0, 0), (1, 1, 0)],
489 [(0, 1, 0), (1, 0, 0)],
490 [(0, 0, 1), (1, 1, 1)],
491 [(0, 1, 1), (1, 0, 1)],
492 ]
493 basis_vectors = [
494 [-0.5, 0.5, 0.5],
495 [ 0.5, -0.5, 0.5],
496 [ 0.5, 0.5, -0.5]
497 ]
498 elif cell[1] == 'M': # Base-centered Monoclinic (底心単斜)
499 a, b, c, alpha, beta, gamma = 5.0, 7.0, 6.0, 90.0, 110.0, 90.0
500 lattice_points = [
501 [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1],
502 [0, 1, 1], [1, 0, 1], [1, 1, 0], [1, 1, 1],
503 [0.5, 0.0, 0.5],
504 [0.5, 1.0, 0.5]
505 ]
506 dashed_lines = [
507 [(0, 0, 0), (1, 0, 1)],
508 [(0, 0, 1), (1, 0, 0)],
509 [(0, 1, 0), (1, 1, 1)],
510 [(0, 1, 1), (1, 1, 0)],
511 ]
512 basis_vectors = [
513 [-0.5, 0.5, 0.5],
514 [ 0.5, -0.5, 0.5],
515 [ 0.5, 0.5, -0.5]
516 ]
517
518 else:
519 print(f"\nError: Invalid cell type [{cell}]\n")
520 return
521
522 draw_unit_cell_with_lattice(
523 a=a, b=b, c=c, alpha=alpha, beta=beta, gamma=gamma,
524 lattice_points=lattice_points,
525 dashed_lines=dashed_lines,
526 basis_vectors=basis_vectors
527 )
528
529
530if __name__ == "__main__":
531 main()