"""
概要: 単位格子と逆格子単位格子を描画するスクリプト。
詳細説明:
本スクリプトは、指定された格子定数と原子サイト情報に基づいて、実空間および逆空間の単位格子と原子配置を3Dで描画します。
格子定数から格子ベクトル、計量テンソル、単位格子の体積を計算し、その結果も表示します。
逆格子についても同様の計算を行い、結果を表示します。
生成された原子サイトは、指定された描画範囲に基づいて重複を避けながらリストに追加されます。
最終的に、matplotlibを用いて実空間の単位格子、逆空間の単位格子、および原子配置を3Dプロットとして可視化し、表示します。
関連リンク: :doc:`crystal_draw_cell_usage`
"""
import sys
import os
import copy
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from tkcrystalbase import *
# Lattice parameters (angstrom and degree)
lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
#lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
sites = [
['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.5, 0.5])]
,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.0, 0.5])]
,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.5, 0.0])]
,['Cl', 'Cl1', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
,['Cl', 'Cl2', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
,['Cl', 'Cl3', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
,['Cl', 'Cl4', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
]
# Coefficient for atomic size to plot
kr = 100.0
# Distance to judge identical atom site, in angstrom
rmin = 0.1
# Coefficient to plot reciprocal unit cell w.r.t. real space unit cell
kRUC = 0.8
# Range of unit cells to draw crystal structure
nrange = [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]]
# Figure configuration
figsize = (12, 12)
[ドキュメント]
def main():
"""
概要: 結晶の単位格子と逆格子単位格子を計算し、3Dプロットで表示する。
詳細説明:
本関数は、グローバル変数 `lattice_parameters` と `sites` で定義された結晶情報を使用します。
1. 格子定数から格子ベクトル、計量テンソル、単位格子の体積を計算し、コンソールに表示します。
2. 逆格子についても同様の計算（逆格子ベクトル、逆格子定数、逆格子計量テンソル、逆格子単位格子の体積）を行い、コンソールに表示します。
3. グローバル変数 `nrange` で指定された範囲に基づいて、結晶内のすべての原子サイトを生成し、重複を避けて `allsites` リストに格納します。
4. matplotlibを用いて、実空間の単位格子、逆空間の単位格子、および原子配置を3Dで描画します。
5. プロットのアスペクト比を調整し、表示します。
引数: なし
戻り値: なし
"""
print("")
print("Lattice parameters:", lattice_parameters)
aij = cal_lattice_vectors(lattice_parameters)
print("Lattice vectors:")
print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
inf = cal_metrics(lattice_parameters)
gij = inf['gij']
print("Metric tensor:")
print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
volume = cal_volume(aij)
print("Volume: {:12.4g} A^3".format(volume))
print("")
print("Unit cell volume: {:12.4g} A^3".format(volume))
Raij = cal_reciprocal_lattice_vectors(aij)
Rlatt = cal_reciprocal_lattice_parameters(Raij)
Rinf = cal_metrics(Rlatt)
Rgij = Rinf['gij']
print("Reciprocal lattice parameters:", Rlatt)
print("Reciprocal lattice vectors:")
print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
print("Reciprocal lattice metric tensor:")
print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
Rvolume = cal_volume(Raij)
print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
allsites = []
for site in sites:
name, label, z, M, q, r, color, pos = copy.deepcopy(site)
pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
if -0.1 <= posn[0] <= 1.1 and -0.1 <= posn[1] <= 1.1 and -0.1 <= posn[2] <= 1.1:
add_site(allsites, [name, label, z, M, q, r, color, posn], gij, rmin)
print("")
print("All sites to draw:")
for s in allsites:
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]))
fig = plt.figure(figsize = figsize)
ax = fig.add_subplot(111, projection='3d')
# Real space unit cell
draw_unitcell(ax, allsites, aij, nrange, kr, linecolor = 'black')
# Reciprocal space unit cell
k = max([*aij[0], *aij[1], *aij[2]]) / max([*Raij[0], *Raij[1], *Raij[2]]) * kRUC
kRaij = np.empty([3, 3])
for i in range(3):
for j in range(3):
kRaij[i][j] = k * Raij[i][j]
draw_unitcell(ax, None, kRaij, nrange, linecolor = 'red')
# Note: set_aspect() is not implemented for 3D plots
# ax.set_aspect('equal','box')
xlim =ax.get_xlim()
ylim =ax.get_ylim()
zlim =ax.get_zlim()
lim = [min([xlim[0], ylim[0], zlim[0]]), max([xlim[1], ylim[1], zlim[1]])]
ax.set_xlim(lim)
ax.set_ylim(lim)
ax.set_zlim(lim)
# ax.set_xticks(np.linspace(*lim, 0))
# ax.set_yticks(np.linspace(*lim, 0))
# ax.set_zticks(np.linspace(*lim, 0))
plt.show()
print("")
exit()
if __name__ == '__main__':
main()