"""
単位格子を異なる基底ベクトルに変換し、変換前後の単位格子を描画するスクリプトです。
このスクリプトは、指定された結晶構造の単位格子を、異なる基底ベクトル（例えばプリミティブセル）を持つ新しい単位格子に変換します。
変換前と変換後の単位格子の格子パラメータ、格子ベクトル、体積、原子サイト情報を計算して表示し、
さらにmatplotlibの3Dプロット機能を使用して両方の単位格子と原子配置を視覚的に表示します。
コマンドライン引数で結晶の種類、変換モード、原子の描画サイズを調整できます。
Requirement: tkcrystalbase.py
:doc:`crystal_convert_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 *
# Crystal name
crystal_name = 'FCC'
# Conversion mode
conversion_mode = 'FCCPrim'
# 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)
crystal_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])]
]
# Coefficient for atomic size to plot
kr = 1000.0
# Distance to judge identical atom site, in angstrom
rmin = 0.1
# 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)
# Atom color for converted cell
converted_atom_color = 'gray'
[ドキュメント]
def get_crystal(name):
"""
指定された結晶名に対応する格子パラメータとサイト情報を返します。
サポートされている結晶タイプはFCC、BCC、Hex、Rhombです。
:param name: str, 結晶名 ('FCC', 'BCC', 'Hex', 'Rhomb')
:returns: tuple, (lattice_parameters, crystal_sites)
- lattice_parameters: list, 格子定数 [a, b, c, alpha, beta, gamma]
- crystal_sites: list, 原子サイト情報リスト
"""
if name == 'FCC':
return lattice_parameters, crystal_sites
if name == 'BCC':
return lattice_parameters, \
[ ['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.5, 0.5, 0.5])] ]
if name == 'Hex':
return [ 5.62, 5.62, 4.00, 90.0, 90.0, 120.0], \
[ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])] ]
if name == 'Rhomb':
return [ 5.62, 5.62, 5.62, 70.0, 70.0, 70.0], \
[ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])] ]
[ドキュメント]
def get_conversion_matrix(key):
"""
指定された変換キーに対応する変換行列を返します。
FCCプリミティブセル、BCCプリミティブセル、面心格子、斜方晶系-菱面体晶、
菱面体晶-六方晶などの変換行列をサポートします。
:param key: str, 変換モードのキー ('FCCPrim', 'BCCPrim', 'ACenterPrim', 'BCenterPrim', 'CCenterPrim', 'RhombHex', 'HexRhomb', 'HexOrtho')
:returns: numpy.ndarray, 3x3の変換行列
"""
if key == 'FCCPrim':
return np.array([ [ 0.5, 0.5, 0],
[ 0, 0.5, 0.5],
[ 0.5, 0, 0.5] ])
elif key == 'BCCPrim':
return np.array([ [-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5] ])
elif key == 'ACenterPrim':
return np.array([ [ 1.0, 0.0, 0.0],
[ 0.0, 0.5, 0.5],
[ 0.0,-0.5, 0.5] ])
elif key == 'BCenterPrim':
return np.array([ [ 0.5, 0.0, 0.5],
[ 0.0, 1.0, 0.0],
[-0.5, 0.0, 0.5] ])
elif key == 'CCenterPrim':
return np.array([ [ 0.5, 0.5, 0.0],
[-0.5, 0.5, 0.0],
[ 0.0, 0.0, 1.0] ])
elif key == 'RhombHex':
return np.array([ [ 1.0, -1.0, 0.0],
[ 0.0, 1.0, -1.0],
[ 1.0, 1.0, 1.0] ])
elif key == 'HexRhomb':
return np.array([ [ 2.0/3, 1.0/3, 1.0/3],
[-1.0/3, 1.0/3, 1.0/3],
[-1.0/3,-2.0/3, 1.0/3] ])
elif key == 'HexOrtho':
return np.array([ [ 1.0, 0.0, 0.0],
[ 1.0, 2.0, 0.0],
[ 0.0, 0.0, 1.0] ])
[ドキュメント]
def usage():
"""
スクリプトの正しい使用方法を標準出力に表示します。
期待されるコマンドライン引数とその説明、および使用例を示します。
:param: None
:returns: None
"""
print("")
print("Usage: python {} crystal_name conversion_mode kRatom".format(argv[0]))
print(" crystal_name : 'FCC', 'BCC', 'Hex', 'Rhomb'")
print(" conversion_mode: 'FCCPrim', 'BCCPrim', 'HexOrtho', 'HexRhomb', 'RhombHex'")
print(" kRatom : Coefficient of atomic radius to draw")
print("ex: python {} {} {} {}".format(argv[0], crystal_name, conversion_mode, kr))
print("")
[ドキュメント]
def terminate():
"""
usageメッセージを表示し、スクリプトを終了します。
エラー発生時や不正な引数が渡された際に呼び出され、ユーザーに正しい使い方を案内した後にプログラムを終了させます。
:param: None
:returns: None
"""
usage()
exit()
[ドキュメント]
def main():
"""
単位格子の変換と描画のメイン処理を実行します。
結晶情報と変換行列を取得し、変換前後の格子パラメータ、格子ベクトル、体積、原子サイト情報を計算して標準出力に表示します。
その後、変換前後の単位格子と原子配置をmatplotlibの3Dプロットで可視化します。
`tkcrystalbase` モジュール内のヘルパー関数を利用して、幾何学的計算と描画を行います。
:param: None
:returns: None
"""
latt, sites = get_crystal(crystal_name)
tij = get_conversion_matrix(conversion_mode)
print("")
print("Crystal name:", crystal_name)
print(" Lattice parameters: {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g}".format(*latt))
aij = cal_lattice_vectors(latt)
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(latt)
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("Sites:")
for s in sites:
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]))
print("")
print("Conversion mode:", conversion_mode)
print("(tij) = |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[0][0], tij[0][1], tij[0][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[1][0], tij[1][1], tij[1][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[2][0], tij[2][1], tij[2][2]))
acij = convert_lattice_vectors(aij, tij)
print(" Converted lattice vectors:")
print(" acx: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[0][0], acij[0][1], acij[0][2]))
print(" acy: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[1][0], acij[1][1], acij[1][2]))
print(" acz: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[2][0], acij[2][1], acij[2][2]))
lattc = cal_lattice_parameters_from_lattice_vectors(acij)
print(" Lattice parameters: {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g}".format(*lattc))
infc = cal_metrics(lattc)
gcij = infc['gij']
print(" Metric tensor:")
print(" gcij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[2]))
volumec = cal_volume(acij)
print(" Volumec: {:12.4g} A^3".format(volumec))
print(" Sites:")
#tij_x2xc = get_conversion_matrix_from_tij(tij, 'x2xc')
tij_x2xc = np.linalg.inv(tij).transpose()
print(" (tij:x to xc) = |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[0][0], tij_x2xc[0][1], tij_x2xc[0][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[1][0], tij_x2xc[1][1], tij_x2xc[1][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[2][0], tij_x2xc[2][1], tij_x2xc[2][2]))
print(" Converted sites:")
csites = []
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):
posc = convert_fractional_coordinates_by_tij([pos01[0] + ix, pos01[1] + iy, pos01[2] + iz], tij)
posn = [reduce01(posc[0]), reduce01(posc[1]), reduce01(posc[2])]
if add_site(csites, [name, label, z, M, q, r * 0.3, color, posn], gcij, rmin):
print(" {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(name, label, *posn, q))
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)
allcsites = []
for site in csites:
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(allcsites, [name, label, z, M, q, r, color, posn], gcij, 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]))
print("")
print("All converted sites to draw:")
for s in allcsites:
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')
draw_unitcell(ax, allcsites, acij, nrange, kr, linecolor = 'blue', atomcolor = converted_atom_color)
# 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("")
terminate()
if __name__ == '__main__':
argv = sys.argv
narg = len(argv)
if narg >= 2:
crystal_name = argv[1]
if narg >= 3:
conversion_mode = argv[2]
if narg >= 4:
kr = float(argv[3])
main()