crystal_draw_cell.py ダウンロード/コピー
crystal_draw_cell.py
crystal_draw_cell.py
1"""
2概要: 単位格子と逆格子単位格子を描画するスクリプト。
3詳細説明:
4 本スクリプトは、指定された格子定数と原子サイト情報に基づいて、実空間および逆空間の単位格子と原子配置を3Dで描画します。
5 格子定数から格子ベクトル、計量テンソル、単位格子の体積を計算し、その結果も表示します。
6 逆格子についても同様の計算を行い、結果を表示します。
7 生成された原子サイトは、指定された描画範囲に基づいて重複を避けながらリストに追加されます。
8 最終的に、matplotlibを用いて実空間の単位格子、逆空間の単位格子、および原子配置を3Dプロットとして可視化し、表示します。
9関連リンク: :doc:`crystal_draw_cell_usage`
10"""
11
12import sys
13import os
14import copy
15from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
16import numpy as np
17from numpy import linalg as la
18
19from tkcrystalbase import *
20
21
22# Lattice parameters (angstrom and degree)
23lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
24#lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
25
26# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
27sites = [
28 ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
29 ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.5, 0.5])]
30 ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.0, 0.5])]
31 ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.5, 0.0])]
32 ,['Cl', 'Cl1', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
33 ,['Cl', 'Cl2', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
34 ,['Cl', 'Cl3', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
35 ,['Cl', 'Cl4', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
36 ]
37
38# Coefficient for atomic size to plot
39kr = 100.0
40
41# Distance to judge identical atom site, in angstrom
42rmin = 0.1
43
44# Coefficient to plot reciprocal unit cell w.r.t. real space unit cell
45kRUC = 0.8
46
47# Range of unit cells to draw crystal structure
48nrange = [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]]
49
50# Figure configuration
51figsize = (12, 12)
52
53
54def main():
55 """
56 概要: 結晶の単位格子と逆格子単位格子を計算し、3Dプロットで表示する。
57 詳細説明:
58 本関数は、グローバル変数 `lattice_parameters` と `sites` で定義された結晶情報を使用します。
59 1. 格子定数から格子ベクトル、計量テンソル、単位格子の体積を計算し、コンソールに表示します。
60 2. 逆格子についても同様の計算（逆格子ベクトル、逆格子定数、逆格子計量テンソル、逆格子単位格子の体積）を行い、コンソールに表示します。
61 3. グローバル変数 `nrange` で指定された範囲に基づいて、結晶内のすべての原子サイトを生成し、重複を避けて `allsites` リストに格納します。
62 4. matplotlibを用いて、実空間の単位格子、逆空間の単位格子、および原子配置を3Dで描画します。
63 5. プロットのアスペクト比を調整し、表示します。
64 引数: なし
65 戻り値: なし
66 """
67 print("")
68 print("Lattice parameters:", lattice_parameters)
69 aij = cal_lattice_vectors(lattice_parameters)
70 print("Lattice vectors:")
71 print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
72 print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
73 print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
74 inf = cal_metrics(lattice_parameters)
75 gij = inf['gij']
76 print("Metric tensor:")
77 print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
78 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
79 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
80 volume = cal_volume(aij)
81 print("Volume: {:12.4g} A^3".format(volume))
82
83 print("")
84 print("Unit cell volume: {:12.4g} A^3".format(volume))
85 Raij = cal_reciprocal_lattice_vectors(aij)
86 Rlatt = cal_reciprocal_lattice_parameters(Raij)
87 Rinf = cal_metrics(Rlatt)
88 Rgij = Rinf['gij']
89 print("Reciprocal lattice parameters:", Rlatt)
90 print("Reciprocal lattice vectors:")
91 print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
92 print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
93 print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
94 print("Reciprocal lattice metric tensor:")
95 print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
96 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
97 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
98 Rvolume = cal_volume(Raij)
99 print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
100
101 allsites = []
102 for site in sites:
103 name, label, z, M, q, r, color, pos = copy.deepcopy(site)
104 pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
105 for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
106 for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
107 for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
108 posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
109 if -0.1 <= posn[0] <= 1.1 and -0.1 <= posn[1] <= 1.1 and -0.1 <= posn[2] <= 1.1:
110 add_site(allsites, [name, label, z, M, q, r, color, posn], gij, rmin)
111
112 print("")
113 print("All sites to draw:")
114 for s in allsites:
115 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]))
116
117
118 fig = plt.figure(figsize = figsize)
119 ax = fig.add_subplot(111, projection='3d')
120
121# Real space unit cell
122 draw_unitcell(ax, allsites, aij, nrange, kr, linecolor = 'black')
123
124# Reciprocal space unit cell
125 k = max([*aij[0], *aij[1], *aij[2]]) / max([*Raij[0], *Raij[1], *Raij[2]]) * kRUC
126 kRaij = np.empty([3, 3])
127 for i in range(3):
128 for j in range(3):
129 kRaij[i][j] = k * Raij[i][j]
130 draw_unitcell(ax, None, kRaij, nrange, linecolor = 'red')
131
132# Note: set_aspect() is not implemented for 3D plots
133# ax.set_aspect('equal','box')
134 xlim =ax.get_xlim()
135 ylim =ax.get_ylim()
136 zlim =ax.get_zlim()
137 lim = [min([xlim[0], ylim[0], zlim[0]]), max([xlim[1], ylim[1], zlim[1]])]
138 ax.set_xlim(lim)
139 ax.set_ylim(lim)
140 ax.set_zlim(lim)
141# ax.set_xticks(np.linspace(*lim, 0))
142# ax.set_yticks(np.linspace(*lim, 0))
143# ax.set_zticks(np.linspace(*lim, 0))
144
145 plt.show()
146
147 print("")
148 exit()
149
150
151if __name__ == '__main__':
152 main()