draw_unit_cells.py ダウンロード/コピー

draw_unit_cells.py をダウンロード

draw_unit_cells.py
draw_unit_cells.py
  1"""
  2単位格子および超格子を描画するためのスクリプト。
  3
  4概要:
  5    結晶学における単位格子（unit cell）や、その繰り返しからなる超格子（supercell）を3Dで可視化します。
  6
  7詳細説明:
  8    格子定数と角度を入力として、3D空間における単位格子を構築し、指定された数だけ繰り返して描画します。
  9    基本ベクトルを計算し、平行六面体を積み重ねることで超格子の形状を表現します。
 10
 11関連リンク: draw_unit_cells_usage
 12"""
 13
 14import numpy as np
 15import matplotlib.pyplot as plt
 16from mpl_toolkits.mplot3d import Axes3D
 17
 18#===================
 19# デフォルト設定
 20#===================
 21ELEV_DEF = 5.90
 22AZIM_DEF = -67.55
 23AXIS_FONT_SIZE_DEF = 24
 24DRAW_LATTICE_VECTORS_DEF = False
 25
 26def lattice_vectors(a, b, c, alpha, beta, gamma):
 27    """
 28    概要:
 29        結晶の格子定数と角度から基本格子ベクトルを計算します。
 30
 31    詳細説明:
 32        直交座標系における3つの基本格子ベクトル va, vb, vc を導出します。
 33        va はX軸に沿い、vb はXY平面にあり、vc は3次元空間に配置されるように定義されます。
 34
 35    引数:
 36        :param a: 格子定数aの長さ。
 37        :type a: float
 38        :param b: 格子定数bの長さ。
 39        :type b: float
 40        :param c: 格子定数cの長さ。
 41        :type c: float
 42        :param alpha: 角度α（度）。vb と vc のなす角。
 43        :type alpha: float
 44        :param beta: 角度β（度）。va と vc のなす角。
 45        :type beta: float
 46        :param gamma: 角度γ（度）。va と vb のなす角。
 47        :type gamma: float
 48    戻り値:
 49        :returns: 3つの基本格子ベクトル va, vb, vc。
 50        :rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]
 51    """
 52    alpha_r, beta_r, gamma_r = np.radians([alpha, beta, gamma])
 53    va = np.array([a, 0, 0])
 54    vb = np.array([b * np.cos(gamma_r), b * np.sin(gamma_r), 0])
 55    cx = c * np.cos(beta_r)
 56    cy = c * (np.cos(alpha_r) - np.cos(beta_r) * np.cos(gamma_r)) / np.sin(gamma_r)
 57    cz = np.sqrt(max(0, c**2 - cx**2 - cy**2))
 58    vc = np.array([cx, cy, cz])
 59    return va, vb, vc
 60
 61def set_equal_aspect(ax, points):
 62    """
 63    概要:
 64        Matplotlib 3Dプロットのアスペクト比を全ての軸で等しく設定します。
 65
 66    引数:
 67        :param ax: 3DプロットのAxesオブジェクト。
 68        :type ax: matplotlib.axes.Axes
 69        :param points: 表示範囲を決定するための点群。
 70        :type points: numpy.ndarray
 71    戻り値:
 72        :returns: None
 73        :rtype: None
 74    """
 75    xlim = [np.min(points[:,0]), np.max(points[:,0])]
 76    ylim = [np.min(points[:,1]), np.max(points[:,1])]
 77    zlim = [np.min(points[:,2]), np.max(points[:,2])]
 78    max_range = max(
 79        xlim[1] - xlim[0],
 80        ylim[1] - ylim[0],
 81        zlim[1] - zlim[0]
 82    ) / 2
 83
 84    mid_x = np.mean(xlim)
 85    mid_y = np.mean(ylim)
 86    mid_z = np.mean(zlim)
 87
 88    ax.set_xlim(mid_x - max_range, mid_x + max_range)
 89    ax.set_ylim(mid_y - max_range, mid_y + max_range)
 90    ax.set_zlim(mid_z - max_range, mid_z + max_range)
 91    ax.set_box_aspect([1,1,1])
 92
 93def draw_cell(ax, va, vb, vc, origin, linewidth=0.5):
 94    """
 95    概要:
 96        指定された基本格子ベクトルと原点に基づいて単一の単位格子を描画します。
 97
 98    引数:
 99        :param ax: 3DプロットのAxesオブジェクト。
100        :type ax: matplotlib.axes.Axes
101        :param va: 第1基本格子ベクトル。
102        :type va: numpy.ndarray
103        :param vb: 第2基本格子ベクトル。
104        :type vb: numpy.ndarray
105        :param vc: 第3基本格子ベクトル。
106        :type vc: numpy.ndarray
107        :param origin: 描画開始点。
108        :type origin: numpy.ndarray
109        :param linewidth: 線の太さ。
110        :type linewidth: float
111    戻り値:
112        :returns: None
113        :rtype: None
114    """
115    p0 = origin
116    p1 = origin + va
117    p2 = origin + vb
118    p3 = origin + vc
119    p4 = origin + va + vb
120    p5 = origin + va + vc
121    p6 = origin + vb + vc
122    p7 = origin + va + vb + vc
123    verts = [p0, p1, p2, p3, p4, p5, p6, p7]
124
125    edges = [
126        (0,1), (0,2), (0,3),
127        (1,4), (1,5),
128        (2,4), (2,6),
129        (3,5), (3,6),
130        (4,7), (5,7), (6,7)
131    ]
132    for i,j in edges:
133        ax.plot(*zip(verts[i], verts[j]), color='black', linewidth=linewidth)
134
135def draw_supercell_plot(a, b, c, alpha, beta, gamma, nx=3, ny=3, nz=3, 
136                        elev=ELEV_DEF, azim=AZIM_DEF, 
137                        draw_lattice_vectors=DRAW_LATTICE_VECTORS_DEF, 
138                        font_size=AXIS_FONT_SIZE_DEF):
139    """
140    概要:
141        指定された格子定数と角度に基づき、超格子を描画し、画像として保存して表示します。
142
143    引数:
144        :param a: 格子定数aの長さ。
145        :type a: float
146        :param b: 格子定数bの長さ。
147        :type b: float
148        :param c: 格子定数cの長さ。
149        :type c: float
150        :param alpha: 角度α（度）。
151        :type alpha: float
152        :param beta: 角度β（度）。
153        :type beta: float
154        :param gamma: 角度γ（度）。
155        :type gamma: float
156        :param nx: X方向の単位格子の繰り返し数。
157        :type nx: int
158        :param ny: Y方向の単位格子の繰り返し数。
159        :type ny: int
160        :param nz: Z方向の単位格子の繰り返し数。
161        :type nz: int
162        :param elev: 3Dプロットの視点仰角。
163        :type elev: float
164        :param azim: 3Dプロットの視点方位角。
165        :type azim: float
166        :param draw_lattice_vectors: 基本格子ベクトルを描画するかどうかのフラグ。
167        :type draw_lattice_vectors: bool
168        :param font_size: 軸ラベルのフォントサイズ。
169        :type font_size: int
170    戻り値:
171        :returns: None
172        :rtype: None
173    """
174    va, vb, vc = lattice_vectors(a, b, c, alpha, beta, gamma)
175    fig = plt.figure(figsize=(8,6))
176    ax = fig.add_subplot(111, projection='3d')
177    ax.view_init(elev=elev, azim=azim)
178
179    # 角度確認用のイベントハンドラ
180    def on_draw(event):
181        curr_elev = ax.elev
182        curr_azim = ax.azim
183        print(f"Current elev: {curr_elev:.2f}, azim: {curr_azim:.2f}")
184    fig.canvas.mpl_connect('draw_event', on_draw)
185
186    # 超格子の全セルの描画（細い黒線）
187    for i in range(nx):
188        for j in range(ny):
189            for k in range(nz):
190                offset = i * va + j * vb + k * vc
191                draw_cell(ax, va, vb, vc, offset, linewidth=0.5)
192
193    # 原点にあるセルの強調（太い黒線）
194    draw_cell(ax, va, vb, vc, np.array([0,0,0]), linewidth=2)
195
196    # 基本格子ベクトルの描画（オプション）
197    if draw_lattice_vectors:
198        for vec, label in zip([va, vb, vc], ['a', 'b', 'c']):
199            ax.quiver(0, 0, 0, *vec, color='blue', linewidth=2, arrow_length_ratio=0.1)
200            ax.text(*vec * 1.05, label, fontsize=font_size, color='blue')
201
202    # 不要な座標軸やグリッドを非表示にする
203    ax.set_xticks([]); ax.set_yticks([]); ax.set_zticks([])
204    ax.set_xticklabels([]); ax.set_yticklabels([]); ax.set_zticklabels([])
205    ax.set_xlabel(''); ax.set_ylabel(''); ax.set_zlabel('')
206    ax.grid(False)
207    ax.xaxis.pane.set_visible(False)
208    ax.yaxis.pane.set_visible(False)
209    ax.zaxis.pane.set_visible(False)
210    ax.xaxis.line.set_color((0,0,0,0))
211    ax.yaxis.line.set_color((0,0,0,0))
212    ax.zaxis.line.set_color((0,0,0,0))
213    ax.set_facecolor((1,1,1,0))
214
215    # 全格子点の計算と描画（小さな●）
216    all_points = []
217    for i in range(nx+1):
218        for j in range(ny+1):
219            for k in range(nz+1):
220                all_points.append(i*va + j*vb + k*vc)
221    all_points = np.array(all_points)
222    ax.scatter(all_points[:,0], all_points[:,1], all_points[:,2], color='black', s=20)
223
224    # アスペクト比の調整
225    set_equal_aspect(ax, all_points)
226
227    plt.tight_layout()
228    plt.savefig("unit_cells.png", dpi=300, bbox_inches='tight', transparent=True)
229    plt.show()
230
231def main():
232    """
233    概要:
234        メインの実行フロー。特定のパラメータで超格子を描画します。
235
236    戻り値:
237        :returns: None
238        :rtype: None
239    """
240    params = {
241        'a': 5.0, 'b': 5.5, 'c': 4.5,
242        'alpha': 80, 'beta': 70, 'gamma': 100,
243        'nx': 3, 'ny': 3, 'nz': 3
244    }
245    
246    print(f"Drawing supercell: {params['nx']}x{params['ny']}x{params['nz']}")
247    draw_supercell_plot(**params)
248
249if __name__ == "__main__":
250    main()