tkcrystalbase.py ダウンロード/コピー
tkcrystalbase.py
tkcrystalbase.py
1"""
2tkcrystalbase.py
3
4概要: 結晶構造計算と描画に関連するユーティリティ関数を提供するモジュール。
5詳細説明:
6 このモジュールは、結晶の格子定数から実空間および逆空間の格子ベクトル、
7 メトリックテンソル、単位格子の体積などを計算する機能を提供します。
8 また、原子の分数座標とデカルト座標間の変換、原子間距離や角度の計算、
9 そしてmatplotlibを使用した3D描画ヘルパー関数を含んでいます。
10 物理定数や一般的な格子パラメータ、サイト情報もここで定義されています。
11
12関連リンク:
13 :doc:`tkcrystalbase_usage`
14"""
15
16from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
17import numpy as np
18from numpy import linalg as la
19from pprint import pprint
20
21pi = 3.14159265358979323846
22pi2 = pi + pi
23torad = 0.01745329251944 # rad/deg";
24todeg = 57.29577951472 # deg/rad";
25basee = 2.71828183
26
27h = 6.6260755e-34 # Js";
28h_bar = 1.05459e-34 # "Js";
29hbar = h_bar
30c = 2.99792458e8 # m/s";
31e = 1.60218e-19 # C";
32me = 9.1093897e-31 # kg";
33mp = 1.6726231e-27 # kg";
34mn = 1.67495e-27 # kg";
35u0 = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
36e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
37e2_4pie0 = 2.30711e-28 # Nm<sup>2</sup>";
38a0 = 5.29177e-11 # m";
39kB = 1.380658e-23 # JK<sup>-1</sup>";
40NA = 6.0221367e23 # mol<sup>-1</sup>";
41R = 8.31451 # J/K/mol";
42F = 96485.3 # C/mol";
43g = 9.81 # m/s2";
44
45
46# Lattice parameters (angstrom and degree)
47#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
48lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
49
50# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
51sites = [
52 ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
53 ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.5, 0.5])]
54 ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.0, 0.5])]
55 ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.5, 0.0])]
56 ,['Cl', 'Cl1', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
57 ,['Cl', 'Cl2', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
58 ,['Cl', 'Cl3', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
59 ,['Cl', 'Cl4', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
60 ]
61
62
63def reduce01(x):
64 """
65 概要: 数値を0以上1未満の範囲に変換する。
66 詳細説明:
67 入力された数値の小数部分を返します。
68 例えば、3.75は0.75に、-2.3は0.7に変換されます。
69 :param x: float: 変換する数値。
70 :returns: float: 0以上1未満の範囲に変換された値。
71 """
72 return x - int(x)
73
74 return x
75
76def round01(x):
77 """
78 概要: 数値を特定の閾値に基づいて0または1に丸めるか、0以上1未満の範囲に変換する。
79 詳細説明:
80 `x` が1.0に近い（絶対差が0.0002未満）場合は1.0を返します。
81 `x` が0.0に近い（絶対差が0.0002未満）場合は0.0を返します。
82 それ以外の場合、`x > 1.0` なら `x - int(x)` を返します。
83 `x < 1.0` なら `x - int(x) + 1.0` を返します。
84 この関数は、負の数や1を超える数値を0以上1未満の範囲に正規化するために使用されることがあります。
85 :param x: float: 丸める、または変換する数値。
86 :returns: float: 丸められた、または0以上1未満の範囲に変換された値。
87 """
88 if abs(x - 1.0) < 0.0002:
89 return 1.0
90 if abs(x) < 0.0002:
91 return 0.0
92 if x > 1.0:
93 return x - int(x)
94 if x < 1.0:
95 return x - int(x) + 1.0
96 return x
97
98def round_parameter(x, tol):
99 """
100 概要: 数値を指定された許容範囲 (tolerance) の整数倍に丸める。
101 詳細説明:
102 入力値 `x` を `tol` の整数倍に丸めます。
103 丸め処理は、`(x + 0.1 * tol) / tol` を整数に変換し、`tol` を掛けることで行われます。
104 これにより、例えば `tol=0.5` で `x=1.2` なら `1.0` に、`x=1.3` なら `1.5` に丸められます。
105 :param x: float: 丸める数値。
106 :param tol: float: 丸めの基準となる許容範囲（単位）。
107 :returns: float: `tol` の整数倍に丸められた値。
108 """
109 val = tol * int( (x+0.1*tol) / tol )
110 return val
111
112
113def cal_lattice_vectors(latt):
114 """
115 概要: 単位格子の格子定数から実空間の格子ベクトルを計算する。
116 詳細説明:
117 格子定数 (a, b, c, α, β, γ) を受け取り、
118 直交座標系における3つの格子ベクトル a_vec, b_vec, c_vec を含む3x3行列を計算します。
119 格子ベクトルは、x軸に沿って a_vec が配置されるように構成されます。
120 角度は度数で指定されますが、内部ではラジアンに変換して計算されます。
121 :param latt: list of float: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
122 長さはオングストローム、角度は度で与えられます。
123 :returns: numpy.ndarray: 3x3 の格子ベクトル行列。
124 `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
125 """
126 cosa = cos(torad * latt[3])
127 cosb = cos(torad * latt[4])
128 cosg = cos(torad * latt[5])
129 sing = sin(torad * latt[5])
130
131 aij = np.empty([3, 3], dtype = float)
132 aij[0][0] = latt[0]
133 aij[0][1] = 0.0;
134 aij[0][2] = 0.0;
135 aij[1][0] = latt[1] * cosg
136 aij[1][1] = latt[1] * sing
137 aij[1][2] = 0.0;
138 aij[2][0] = latt[2] * cosb
139 aij[2][1] = latt[2] * (cosa - cosb * cosg) / sing
140 if abs(aij[2][1]) < 1.0e-8:
141 aij[2][1] = 0.0
142 else:
143 aij[2][1] = aij[2][1] / sing
144 aij[2][2] = sqrt(latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1])
145
146 return aij
147
148def cal_metrics(latt):
149 """
150 概要: 格子定数から格子ベクトルとメトリックテンソルを計算する。
151 詳細説明:
152 単位格子の格子定数から実空間の格子ベクトル `aij` を計算し、
153 それらの内積からメトリックテンソル `gij` を構築します。
154 結果は辞書形式で返されます。
155 :param latt: list of float: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
156 長さはオングストローム、角度は度で与えられます。
157 :returns: dict: 計算された格子ベクトル行列 (`'aij'`) とメトリックテンソル (`'gij'`)
158 をキーに持つ辞書。
159 `{'aij': numpy.ndarray, 'gij': numpy.ndarray}` の形式。
160 """
161 inf = {}
162
163 cosa = cos(torad * latt[3])
164 cosb = cos(torad * latt[4])
165 cosg = cos(torad * latt[5])
166
167 aij = cal_lattice_vectors(latt)
168 inf['aij'] = aij
169
170 gij = np.empty([3, 3], dtype = float)
171 for i in range(3):
172 for j in range(i, 3):
173 gij[i][j] = np.dot(aij[i], aij[j])
174 gij[j][i] = gij[i][j]
175 inf['gij'] = gij
176
177 return inf
178
179def cal_volume(aij):
180 """
181 概要: 格子ベクトルから単位格子の体積を計算する。
182 詳細説明:
183 3つの格子ベクトル (a_vec, b_vec, c_vec) のスカラー三重積 (a_vec・(b_vec × c_vec)) を用いて、
184 単位格子の体積を計算します。
185 :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
186 `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
187 :returns: float: 単位格子の体積。
188 """
189 axb = np.cross(aij[0], aij[1]) # Outner product
190 vol = np.dot(axb, aij[2]) # Inner product
191 return vol
192
193def cal_reciprocal_lattice_vectors(aij):
194 """
195 概要: 実空間格子ベクトルから逆空間格子ベクトルを計算する。
196 詳細説明:
197 実空間の格子ベクトル `aij` (a_vec, b_vec, c_vec) を用いて、
198 逆空間の格子ベクトル (a*_vec, b*_vec, c*_vec) を計算します。
199 逆空間ベクトルは、単位格子の体積 V と実空間ベクトルの外積によって定義されます。
200 :param aij: numpy.ndarray: 3x3 の実空間格子ベクトル行列。
201 `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
202 :returns: list of numpy.ndarray: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
203 """
204 V = cal_volume(aij)
205 Ra = np.cross(aij[1], aij[2]) / V
206 Rb = np.cross(aij[2], aij[0]) / V
207 Rc = np.cross(aij[0], aij[1]) / V
208
209 return [Ra, Rb, Rc]
210
211def cal_reciprocal_lattice_parameters(Raij):
212 """
213 概要: 逆空間格子ベクトルから逆空間の格子定数を計算する。
214 詳細説明:
215 逆空間の格子ベクトル `Raij` (a*_vec, b*_vec, c*_vec) を用いて、
216 逆空間の格子定数 (a*, b*, c*, α*, β*, γ*) を計算します。
217 長さはベクトルのノルムから、角度はベクトルの内積から導出されます。
218 :param Raij: list of numpy.ndarray: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
219 :returns: list of float: 逆空間の格子定数 `[Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]` のリスト。
220 """
221 Ra = la.norm(Raij[0])
222 Rb = la.norm(Raij[1])
223 Rc = la.norm(Raij[2])
224 Ralpha = todeg * arccos(np.dot(Raij[1], Raij[2]) / Rb / Rc)
225 Rbeta = todeg * arccos(np.dot(Raij[2], Raij[0]) / Rc / Ra)
226 Rgamma = todeg * arccos(np.dot(Raij[0], Raij[1]) / Ra / Rb)
227
228 return [Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]
229
230def fractional_to_cartesian(pos, aij):
231 """
232 概要: 分数座標からデカルト座標に変換する。
233 詳細説明:
234 単位格子内の原子の分数座標 (x_frac, y_frac, z_frac) と格子ベクトル `aij` を用いて、
235 原子のデカルト座標 (x_cart, y_cart, z_cart) を計算します。
236 :param pos: numpy.ndarray or list of float: 変換する分数座標 `[x_frac, y_frac, z_frac]`。
237 :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
238 `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
239 :returns: tuple of float: 計算されたデカルト座標 `(x_cart, y_cart, z_cart)`。
240 """
241 x = pos[0] * aij[0][0] + pos[1] * aij[1][0] + pos[2] * aij[2][0]
242 y = pos[0] * aij[0][1] + pos[1] * aij[1][1] + pos[2] * aij[2][1]
243 z = pos[0] * aij[0][2] + pos[1] * aij[1][2] + pos[2] * aij[2][2]
244
245 return x, y, z
246
247def distance2(pos0, pos1, gij):
248 """
249 概要: 2つの分数座標間の距離の2乗を計算する。
250 詳細説明:
251 2つの原子の分数座標 `pos0` と `pos1` の間の距離の2乗を、
252 メトリックテンソル `gij` を使用して計算します。
253 この関数は周期境界条件を考慮しません。
254 :param pos0: numpy.ndarray or list of float: 1点目の分数座標 `[x0, y0, z0]`。
255 :param pos1: numpy.ndarray or list of float: 2点目の分数座標 `[x1, y1, z1]`。
256 :param gij: numpy.ndarray: 3x3 のメトリックテンソル。
257 :returns: float: 2点間の距離の2乗。
258 """
259 dx = pos1 - pos0
260# dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
261 r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
262 + 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
263
264 return r2
265
266def distance(pos0, pos1, gij):
267 """
268 概要: 2つの分数座標間の距離を計算する。
269 詳細説明:
270 2つの原子の分数座標 `pos0` と `pos1` の間の距離を、
271 メトリックテンソル `gij` を使用して計算します。
272 この関数は周期境界条件を考慮しません。
273 :param pos0: numpy.ndarray or list of float: 1点目の分数座標 `[x0, y0, z0]`。
274 :param pos1: numpy.ndarray or list of float: 2点目の分数座標 `[x1, y1, z1]`。
275 :param gij: numpy.ndarray: 3x3 のメトリックテンソル。
276 :returns: float: 2点間の距離。
277 """
278 dx = pos1 - pos0
279# dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
280 r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
281 + 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
282 r = sqrt(r2)
283
284 return r
285
286def angle(pos0, pos1, pos2, gij):
287 """
288 概要: 3つの分数座標で定義される角度を計算する。
289 詳細説明:
290 `pos0` を頂点とし、`pos0-pos1` と `pos0-pos2` の2つのベクトルがなす角度を計算します。
291 角度の計算にはメトリックテンソル `gij` が使用されます。
292 結果は度数で返され、180度を超える場合は `360.0 - angle` として調整されます。
293 :param pos0: numpy.ndarray or list of float: 角度の頂点となる分数座標 `[x0, y0, z0]`。
294 :param pos1: numpy.ndarray or list of float: 1つ目のベクトルを定義する分数座標 `[x1, y1, z1]`。
295 :param pos2: numpy.ndarray or list of float: 2つ目のベクトルを定義する分数座標 `[x2, y2, z2]`。
296 :param gij: numpy.ndarray: 3x3 のメトリックテンソル。
297 :returns: float: 3点間で形成される角度（度）。
298 いずれかの距離が0の場合は0.0を返します。
299 """
300 dis01 = distance(pos0, pos1, gij)
301 if dis01 == 0.0:
302 return 0.0
303 dis02 = distance(pos0, pos2, gij)
304 if dis02 == 0.0:
305 return 0.0
306
307 dx01 = pos1 - pos0
308 dx02 = pos2 - pos0
309# dx01 = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
310# dx02 = [pos2[0] - pos0[0], pos2[1] - pos0[1], pos2[2] - pos0[2]]
311 ip = gij[0][0] * dx01[0]*dx02[0] + gij[1][1] * dx01[1]*dx02[1] + gij[2][2] * dx01[2]*dx02[2] \
312 + 2.0 * (gij[0][1] * dx01[0]*dx02[1] + gij[0][2] * dx01[0]*dx02[2] + gij[1][2] * dx01[1]*dx02[2])
313
314 cosa = ip / dis01 / dis02
315 angle = todeg * arccos(cosa)
316 if angle > 180.0:
317 angle = 360.0 - angle
318
319 return angle;
320
321def configure_axis_structure(ax, xrange, yrange, zrange, fontsize = 12, legend_fontsize = 12):
322 """
323 概要: Matplotlibの3D軸オブジェクトの表示設定を行う。
324 詳細説明:
325 MatplotlibのAxes3Dオブジェクト (`ax`) に対して、軸ラベルのフォントサイズ、
326 軸の線、目盛り、グリッド、ペインの表示設定をカスタマイズします。
327 デフォルトでは、軸の線、目盛り、ラベルは非表示に設定され、
328 箱型のアスペクト比で視覚的に整えられます。
329 :param ax: matplotlib.axes.Axes: 設定を行うMatplotlibの3D軸オブジェクト。
330 :param xrange: tuple or list of float: x軸の表示範囲 `(xmin, xmax)`。
331 :param yrange: tuple or list of float: y軸の表示範囲 `(ymin, ymax)`。
332 :param zrange: tuple or list of float: z軸の表示範囲 `(zmin, zmax)`。
333 :param fontsize: int, optional: 軸ラベルおよび目盛りのフォントサイズ。
334 0以下の場合は軸ラベルと目盛りを非表示にします。デフォルトは12。
335 :param legend_fontsize: int, optional: 凡例のフォントサイズ（この関数では直接使用されません）。デフォルトは12。
336 :returns: None
337 """
338 ax.tick_params(axis = 'both', which = 'major', labelsize = fontsize)
339 if fontsize > 0:
340 ax.set_xlabel(f'$x$', fontsize = fontsize, labelpad = -5)
341 ax.set_ylabel(f'$y$', fontsize = fontsize, labelpad = -5)
342 ax.set_zlabel(f'$z$', fontsize = fontsize, labelpad = 0)
343
344# x軸、y軸、z軸の線を非表示にする
345 ax.xaxis.line.set_visible(False)
346 ax.yaxis.line.set_visible(False)
347 ax.zaxis.line.set_visible(False)
348
349 ax.grid(False)
350 ax.xaxis.pane.fill = False
351 ax.yaxis.pane.fill = False
352 ax.zaxis.pane.fill = False
353 ax.xaxis.pane.set_edgecolor('w')
354 ax.yaxis.pane.set_edgecolor('w')
355 ax.zaxis.pane.set_edgecolor('w')
356 ax.xaxis.pane.set_alpha(0)
357 ax.yaxis.pane.set_alpha(0)
358 ax.zaxis.pane.set_alpha(0)
359# ax.w_xaxis.line.set_color('none')
360# ax.w_yaxis.line.set_color('none')
361# ax.w_zaxis.line.set_color('none')
362 ax.set_xticks([])
363 ax.set_yticks([])
364 ax.set_zticks([])
365 ax.set_xticklabels([])
366 ax.set_yticklabels([])
367 ax.set_zticklabels([])
368
369 ax.set_xlim(xrange)
370 ax.set_ylim(yrange)
371 ax.set_zlim(zrange)
372 ax.set_aspect('equal','box')
373# ax.set_box_aspect([xrange[1] - xrange[0], yrange[1] - yrange[0], zrange[1] - zrange[0]])
374
375def draw_box(ax, aij, nrange, color = 'black'):
376 """
377 概要: 3D軸上に単位格子のボックスを描画する。
378 詳細説明:
379 格子ベクトル `aij` で定義される単位格子の12本の辺を、Matplotlibの3D軸 `ax` に描画します。
380 この関数は、単一の単位格子を描画するために設計されており、`nrange` は現在の実装では使用されません。
381 :param ax: matplotlib.axes.Axes: 描画対象のMatplotlib 3D軸オブジェクト。
382 :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
383 `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
384 :param nrange: list of list of int: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
385 この関数では使用されません。
386 :param color: str, optional: ボックスの線の色。デフォルトは 'black'。
387 :returns: None
388 """
389# (0,0,0) -> ax
390 ax.plot([0.0, aij[0][0]],
391 [0.0, aij[0][1]],
392 [0.0, aij[0][2]], color = color)
393# (0,0,0) -> ay
394 ax.plot([0.0, aij[1][0]],
395 [0.0, aij[1][1]],
396 [0.0, aij[1][2]], color = color)
397# (0,0,0) -> az
398 ax.plot([0.0, aij[2][0]],
399 [0.0, aij[2][1]],
400 [0.0, aij[2][2]], color = color)
401
402# ax -> ax + ay
403 ax.plot([aij[0][0], aij[0][0] + aij[1][0]],
404 [aij[0][1], aij[0][1] + aij[1][1]],
405 [aij[0][2], aij[0][2] + aij[1][2]], color = color)
406# ax -> ax + az
407 ax.plot([aij[0][0], aij[0][0] + aij[2][0]],
408 [aij[0][1], aij[0][1] + aij[2][1]],
409 [aij[0][2], aij[0][2] + aij[2][2]], color = color)
410
411# ay -> ay + ax
412 ax.plot([aij[1][0], aij[1][0] + aij[0][0]],
413 [aij[1][1], aij[1][1] + aij[0][1]],
414 [aij[1][2], aij[1][2] + aij[0][2]], color = color)
415# ay -> ay + az
416 ax.plot([aij[1][0], aij[1][0] + aij[2][0]],
417 [aij[1][1], aij[1][1] + aij[2][1]],
418 [aij[1][2], aij[1][2] + aij[2][2]], color = color)
419
420# az -> az + ax
421 ax.plot([aij[2][0], aij[2][0] + aij[0][0]],
422 [aij[2][1], aij[2][1] + aij[0][1]],
423 [aij[2][2], aij[2][2] + aij[0][2]], color = color)
424# az -> ax + ay
425 ax.plot([aij[2][0], aij[2][0] + aij[1][0]],
426 [aij[2][1], aij[2][1] + aij[1][1]],
427 [aij[2][2], aij[2][2] + aij[1][2]], color = color)
428
429# ax + ay -> ax + ay + az
430 ax.plot([aij[0][0] + aij[1][0], aij[0][0] + aij[1][0] + aij[2][0]],
431 [aij[0][1] + aij[1][1], aij[0][1] + aij[1][1] + aij[2][1]],
432 [aij[0][2] + aij[1][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
433
434# ax + az -> ax + ay + az
435 ax.plot([aij[0][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
436 [aij[0][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
437 [aij[0][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
438
439# ay + az -> ax + ay + az
440 ax.plot([aij[1][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
441 [aij[1][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
442 [aij[1][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
443
444def draw_unitcell(ax, sites, aij, nrange, color = 'black', facecolor = 'black', edgecolor = 'white', alpha = 0.7, kr = 1.0):
445 """
446 概要: 3D軸上に単位格子とその内部のサイト原子を描画する。
447 詳細説明:
448 まず `draw_box` 関数を呼び出して単位格子の境界線を描画します。
449 次に、`sites` リストに含まれる各原子について、その分数座標をデカルト座標に変換し、
450 指定された `nrange` の範囲内の周期的な位置に原子を散布図として描画します。
451 原子の描画スタイル（色、透明度、サイズ）は引数でカスタマイズ可能です。
452 :param ax: matplotlib.axes.Axes: 描画対象のMatplotlib 3D軸オブジェクト。
453 :param sites: list of list: サイト情報（原子名、ラベル、原子番号、原子質量、電荷、半径、色、位置）のリスト。
454 各要素は `[name, label, z, M, q, r, color, pos]` の形式。
455 `None` の場合、原子は描画されません。
456 :param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
457 `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
458 :param nrange: list of list of int: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
459 原子の繰り返し構造を描画するために使用されます。
460 :param color: str, optional: 単位格子ボックスの線の色。デフォルトは 'black'。
461 :param facecolor: str, optional: サイト原子の面の描画色。デフォルトは 'black'。
462 :param edgecolor: str, optional: サイト原子のエッジの描画色。デフォルトは 'white'。
463 :param alpha: float, optional: サイト原子の透明度（0.0から1.0）。デフォルトは 0.7。
464 :param kr: float, optional: サイト原子の描画サイズに対するスケール因子。デフォルトは 1.0。
465 :returns: None
466 """
467 draw_box(ax, aij, nrange, color)
468
469 if sites is None: return
470
471 for isite, site in enumerate(sites):
472 name, label, z, M, q, r, color, pos = site
473 pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
474 for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
475 for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
476 for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
477 posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
478 if posn[0] < nrange[0][0] or nrange[0][1] < posn[0] \
479 or posn[1] < nrange[1][0] or nrange[1][1] < posn[1] \
480 or posn[2] < nrange[2][0] or nrange[2][1] < posn[2]:
481 continue
482
483 x, y, z = fractional_to_cartesian(posn, aij)
484 ax.scatter([x], [y], [z], marker = 'o', c = facecolor, edgecolors = edgecolor, alpha = alpha, s = kr *r)
485
486
487def main():
488 """
489 概要: モジュールの主要な計算と結果の表示を実行する。
490 詳細説明:
491 定義された `lattice_parameters` を使用して、実空間および逆空間の
492 格子ベクトル、メトリックテンソル、単位格子の体積を計算します。
493 さらに、サンプルとして原子間距離と角度の計算も行い、
494 その結果を標準出力に表示します。
495 最後にプログラムを終了します。
496 :returns: None
497 """
498 print("")
499 print("Lattice parameters:", lattice_parameters)
500 aij = cal_lattice_vectors(lattice_parameters)
501 print("Lattice vectors:")
502 print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
503 print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
504 print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
505 inf = cal_metrics(lattice_parameters)
506 gij = inf['gij']
507 print("Metric tensor:")
508 print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
509 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
510 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
511 print("")
512 volume = cal_volume(aij)
513 print("Volume: {:12.4g} A^3".format(volume))
514
515 print("")
516 print("Unit cell volume: {:12.4g} A^3".format(volume))
517 Raij = cal_reciprocal_lattice_vectors(aij)
518 Rlatt = cal_reciprocal_lattice_parameters(Raij)
519 Rinf = cal_metrics(Rlatt)
520 Rgij = Rinf['gij']
521 print("Reciprocal lattice parameters:", Rlatt)
522 print("Reciprocal lattice vectors:")
523 print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
524 print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
525 print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
526 print("Reciprocal lattice metric tensor:")
527 print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
528 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
529 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
530 Rvolume = cal_volume(Raij)
531 print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
532
533 print("")
534 print("dis=", distance(np.array([0,0,0]), np.array([1,1,1]), gij))
535 print("angle=", angle (np.array([0,0,0]), np.array([1,1,1]), np.array([1,0,0]), gij))
536 print("angle=", angle (np.array([0,0,0]), np.array([1,0,0]), np.array([0,1,0]), gij))
537
538 print("")
539 exit()
540
541
542if __name__ == '__main__':
543 main()