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

crystal_convert_cell.py をダウンロード

crystal_convert_cell.py
crystal_convert_cell.py
  1"""
  2単位格子を異なる基底ベクトルに変換し、変換前後の単位格子を描画するスクリプトです。
  3
  4概要:
  5    このスクリプトは、指定された結晶構造の単位格子を、異なる基底ベクトル（例えばプリミティブセル）を持つ
  6    新しい単位格子に変換し、その結果を可視化します。
  7詳細説明:
  8    変換前と変換後の単位格子の格子パラメータ、格子ベクトル、体積、原子サイト情報を計算して標準出力に表示します。
  9    さらに、matplotlibの3Dプロット機能を使用して両方の単位格子と原子配置を視覚的に表示します。
 10    コマンドライン引数で結晶の種類、変換モード、原子の描画サイズを調整できます。
 11    このスクリプトは tkcrystalbase.py モジュールに依存しています。
 12関連リンク:
 13    crystal_convert_cell_usage
 14"""
 15import sys
 16import os
 17import copy
 18from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
 19import numpy as np
 20from numpy import linalg as la
 21
 22from tkcrystalbase import *
 23
 24
 25# Crystal name
 26crystal_name = 'FCC'
 27# Conversion mode
 28conversion_mode = 'FCCPrim'
 29
 30# Lattice parameters (angstrom and degree)
 31#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
 32lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
 33
 34# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
 35crystal_sites = [
 36         ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
 37        ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.5, 0.5])]
 38        ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.0, 0.5])]
 39        ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.5, 0.0])]
 40        ]
 41
 42# Coefficient for atomic size to plot
 43kr = 1000.0
 44
 45# Distance to judge identical atom site, in angstrom
 46rmin = 0.1
 47
 48# Range of unit cells to draw crystal structure
 49nrange = [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]]
 50
 51# Figure configuration
 52figsize = (12, 12)
 53
 54# Atom color for converted cell
 55converted_atom_color = 'gray'
 56
 57def get_crystal(name):
 58    """
 59    概要:
 60        指定された結晶名に対応する格子パラメータとサイト情報を返します。
 61    詳細説明:
 62        サポートされている結晶タイプはFCC、BCC、Hex、Rhombです。
 63    引数:
 64        :param name: 結晶名。 'FCC', 'BCC', 'Hex', 'Rhomb' のいずれか。
 65        :type name: str
 66    戻り値:
 67        :returns: 格子定数 (lattice_parameters) と原子サイト情報 (crystal_sites) のタプル。
 68        :rtype: tuple
 69    """
 70    if name == 'FCC':
 71        return lattice_parameters, crystal_sites
 72    if name == 'BCC':
 73        return lattice_parameters, \
 74                [ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
 75                 ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.5, 0.5])] ]
 76    if name == 'Hex':
 77        return [ 5.62, 5.62, 4.00, 90.0, 90.0, 120.0], \
 78                [ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])] ]
 79    if name == 'Rhomb':
 80        return [ 5.62, 5.62, 5.62, 70.0, 70.0, 70.0], \
 81                [ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])] ]
 82
 83def get_conversion_matrix(key):
 84    """
 85    概要:
 86        指定された変換キーに対応する変換行列を返します。
 87    詳細説明:
 88        FCCプリミティブセル、BCCプリミティブセル、面心格子、斜方晶系-菱面体晶、
 89        菱面体晶-六方晶などの様々な結晶学的変換行列をサポートします。
 90    引数:
 91        :param key: 変換モードのキー。 'FCCPrim', 'BCCPrim', 'ACenterPrim', 'BCenterPrim', 'CCenterPrim', 'RhombHex', 'HexRhomb', 'HexOrtho' のいずれか。
 92        :type key: str
 93    戻り値:
 94        :returns: 3x3の変換行列。
 95        :rtype: numpy.ndarray
 96    """
 97    if key == 'FCCPrim':
 98        return np.array([ [ 0.5, 0.5,   0], 
 99                 [ 0,   0.5, 0.5],
100                 [ 0.5,   0, 0.5] ])
101    elif key == 'BCCPrim':
102        return np.array([ [-0.5, 0.5, 0.5], 
103                 [ 0.5,-0.5, 0.5],
104                 [ 0.5, 0.5,-0.5] ])
105    elif key == 'ACenterPrim':
106        return np.array([ [ 1.0, 0.0, 0.0], 
107                 [ 0.0, 0.5, 0.5],
108                 [ 0.0,-0.5, 0.5] ])
109    elif key == 'BCenterPrim':
110        return np.array([ [ 0.5, 0.0, 0.5], 
111                 [ 0.0, 1.0, 0.0],
112                 [-0.5, 0.0, 0.5] ])
113    elif key == 'CCenterPrim':
114        return np.array([ [ 0.5, 0.5, 0.0], 
115                 [-0.5, 0.5, 0.0],
116                 [ 0.0, 0.0, 1.0] ])
117    elif key == 'RhombHex':
118        return np.array([ [ 1.0, -1.0,  0.0], 
119                 [ 0.0,  1.0, -1.0],
120                 [ 1.0,  1.0,  1.0] ])
121    elif key == 'HexRhomb':
122        return np.array([ [ 2.0/3, 1.0/3, 1.0/3],
123                 [-1.0/3, 1.0/3, 1.0/3],
124                 [-1.0/3,-2.0/3, 1.0/3] ])
125    elif key == 'HexOrtho':
126        return np.array([ [ 1.0, 0.0, 0.0], 
127                 [ 1.0, 2.0, 0.0],
128                 [ 0.0, 0.0, 1.0] ])
129
130
131def usage():
132    """
133    概要:
134        スクリプトの正しい使用方法を標準出力に表示します。
135    詳細説明:
136        この関数は、期待されるコマンドライン引数とその説明、および具体的な使用例をユーザーに提示します。
137    戻り値:
138        :returns: なし
139        :rtype: None
140    """
141    print("")
142    print("Usage: python {} crystal_name conversion_mode kRatom".format(argv[0]))
143    print("   crystal_name   : 'FCC',     'BCC',     'Hex',                  'Rhomb'")
144    print("   conversion_mode: 'FCCPrim', 'BCCPrim', 'HexOrtho', 'HexRhomb', 'RhombHex'")
145    print("   kRatom         : Coefficient of atomic radius to draw")
146    print("ex: python {} {} {} {}".format(argv[0], crystal_name, conversion_mode, kr))
147    print("")
148
149def terminate():
150    """
151    概要:
152        usageメッセージを表示し、スクリプトを終了します。
153    詳細説明:
154        エラー発生時や不正な引数が渡された際に呼び出され、ユーザーに正しい使い方を案内した後にプログラムを終了させます。
155    戻り値:
156        :returns: なし
157        :rtype: None
158    """
159    usage()
160    exit()
161
162
163def main():
164    """
165    概要:
166        単位格子の変換と描画のメイン処理を実行します。
167    詳細説明:
168        結晶情報と変換行列を取得し、変換前後の格子パラメータ、格子ベクトル、体積、原子サイト情報を計算して標準出力に表示します。
169        その後、変換前後の単位格子と原子配置をmatplotlibの3Dプロットで可視化します。
170        tkcrystalbase モジュール内のヘルパー関数を利用して、幾何学的計算と描画を行います。
171    戻り値:
172        :returns: なし
173        :rtype: None
174    """
175    latt, sites = get_crystal(crystal_name)
176    tij = get_conversion_matrix(conversion_mode)
177
178    print("")
179    print("Crystal name:", crystal_name)
180    print("  Lattice parameters: {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g}".format(*latt))
181    aij = cal_lattice_vectors(latt)
182    print("Lattice vectors:")
183    print("  ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
184    print("  ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
185    print("  az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
186    inf = cal_metrics(latt)
187    gij = inf['gij']
188    print("Metric tensor:")
189    print("  gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
190    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
191    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
192    volume = cal_volume(aij)
193    print("Volume: {:12.4g} A^3".format(volume))
194    print("Sites:")
195    for s in sites:
196        print("  {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(s[0], s[1], s[7][0], s[7][1], s[7][2], s[4]))
197
198    print("")
199    print("Conversion mode:", conversion_mode)
200    print("(tij) = |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[0][0], tij[0][1], tij[0][2]))
201    print("        |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[1][0], tij[1][1], tij[1][2]))
202    print("        |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[2][0], tij[2][1], tij[2][2]))
203    acij = convert_lattice_vectors(aij, tij)
204    print("  Converted lattice vectors:")
205    print("    acx: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[0][0], acij[0][1], acij[0][2]))
206    print("    acy: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[1][0], acij[1][1], acij[1][2]))
207    print("    acz: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[2][0], acij[2][1], acij[2][2]))
208    lattc = cal_lattice_parameters_from_lattice_vectors(acij)
209    print("    Lattice parameters: {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g}".format(*lattc))
210    infc = cal_metrics(lattc)
211    gcij = infc['gij']
212    print("  Metric tensor:")
213    print("    gcij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[0]))
214    print("          ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[1]))
215    print("          ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[2]))
216    volumec = cal_volume(acij)
217    print("  Volumec: {:12.4g} A^3".format(volumec))
218    print("  Sites:")
219    #tij_x2xc = get_conversion_matrix_from_tij(tij, 'x2xc')
220    tij_x2xc = np.linalg.inv(tij).transpose()
221    print("  (tij:x to xc) = |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[0][0], tij_x2xc[0][1], tij_x2xc[0][2]))
222    print("                  |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[1][0], tij_x2xc[1][1], tij_x2xc[1][2]))
223    print("                  |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[2][0], tij_x2xc[2][1], tij_x2xc[2][2]))
224    print("  Converted sites:")
225    csites = []
226    for site in sites:
227        name, label, z, M, q, r, color, pos = copy.deepcopy(site)
228        pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
229        for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
230         for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
231          for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
232            posc = convert_fractional_coordinates_by_tij([pos01[0] + ix, pos01[1] + iy, pos01[2] + iz], tij)
233            posn = [reduce01(posc[0]), reduce01(posc[1]), reduce01(posc[2])]
234            if add_site(csites, [name, label, z, M, q, r * 0.3, color, posn], gcij, rmin):
235                print("  {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(name, label, *posn, q))
236    
237    allsites = []
238    for site in sites:
239        name, label, z, M, q, r, color, pos = copy.deepcopy(site)
240        pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
241        for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
242         for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
243          for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
244            posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
245            if -0.1 <= posn[0] <= 1.1 and -0.1 <= posn[1] <= 1.1 and -0.1 <= posn[2] <= 1.1:
246                add_site(allsites, [name, label, z, M, q, r, color, posn], gij, rmin)
247
248    allcsites = []
249    for site in csites:
250        name, label, z, M, q, r, color, pos = copy.deepcopy(site)
251        pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
252        for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
253         for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
254          for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
255            posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
256            if -0.1 <= posn[0] <= 1.1 and -0.1 <= posn[1] <= 1.1 and -0.1 <= posn[2] <= 1.1:
257                add_site(allcsites, [name, label, z, M, q, r, color, posn], gcij, rmin)
258
259    print("")
260    print("All sites to draw:")
261    for s in allsites:
262        print("  {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(s[0], s[1], s[7][0], s[7][1], s[7][2], s[4]))
263
264    print("")
265    print("All converted sites to draw:")
266    for s in allcsites:
267        print("  {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(s[0], s[1], s[7][0], s[7][1], s[7][2], s[4]))
268
269
270    fig = plt.figure(figsize = figsize)
271    ax = fig.add_subplot(111, projection='3d')
272
273# Real space unit cell
274    draw_unitcell(ax, allsites,  aij,  nrange, kr, linecolor = 'black')
275    draw_unitcell(ax, allcsites, acij, nrange, kr, linecolor = 'blue', atomcolor = converted_atom_color)
276
277# Note: set_aspect() is not implemented for 3D plots
278#    ax.set_aspect('equal','box')
279    xlim =ax.get_xlim()
280    ylim =ax.get_ylim()
281    zlim =ax.get_zlim()
282    lim = [min([xlim[0], ylim[0], zlim[0]]), max([xlim[1], ylim[1], zlim[1]])]
283    ax.set_xlim(lim)
284    ax.set_ylim(lim)
285    ax.set_zlim(lim)
286#    ax.set_xticks(np.linspace(*lim, 0))
287#    ax.set_yticks(np.linspace(*lim, 0))
288#    ax.set_zticks(np.linspace(*lim, 0))
289
290    plt.show()
291
292    print("")
293    
294
295    terminate()
296
297
298if __name__ == '__main__':
299    argv = sys.argv
300    narg = len(argv)
301    if narg >= 2:
302        crystal_name    = argv[1]
303    if narg >= 3:
304        conversion_mode = argv[2]
305    if narg >= 4:
306        kr = float(argv[3])
307
308    main()