"""
概要: 結晶構造と逆格子構造の単位胞を3Dで描画するスクリプトです。
詳細説明:
このスクリプトは、指定された格子パラメータとサイト情報に基づいて、
実空間における結晶単位胞と、それに対応する逆空間の単位胞を
matplotlibの3Dプロットに描画します。
格子ベクトル、メトリックテンソル、体積などの情報も計算して表示します。
関連リンク: :doc:`crystal_draw_cell_usage`
"""
import sys
import os
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
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
# 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 draw_box(ax, aij, nrange, color = 'black'):
"""
概要: 3Dプロットに結晶単位胞の枠（平行六面体）を描画します。
詳細説明:
3つの格子ベクトルに基づいて、そのベクトルが張る箱の12辺を描画します。
この関数は`nrange`引数を受け取りますが、内部では使用していません。
:param ax: matplotlib.axes.Axesオブジェクト。3Dプロットを描画する軸。
:type ax: mpl_toolkits.mplot3d.axes3d.Axes3D
:param aij: 3x3の行列で、各行が結晶格子ベクトル (ax, ay, az) を表します。
:type aij: numpy.ndarray
:param nrange: 描画範囲を定義するリスト (例: [[min_x, max_x], [min_y, max_y], [min_z, max_z]])。
:type nrange: list
:param color: 描画する辺の色。デフォルトは'black'です。
:type color: str, optional
:returns: なし
:rtype: None
"""
# (0,0,0) -> ax
ax.plot([0.0, aij[0][0]],
[0.0, aij[0][1]],
[0.0, aij[0][2]], color = color)
# (0,0,0) -> ay
ax.plot([0.0, aij[1][0]],
[0.0, aij[1][1]],
[0.0, aij[1][2]], color = color)
# (0,0,0) -> az
ax.plot([0.0, aij[2][0]],
[0.0, aij[2][1]],
[0.0, aij[2][2]], color = color)
# ax -> ax + ay
ax.plot([aij[0][0], aij[0][0] + aij[1][0]],
[aij[0][1], aij[0][1] + aij[1][1]],
[aij[0][2], aij[0][2] + aij[1][2]], color = color)
# ax -> ax + az
ax.plot([aij[0][0], aij[0][0] + aij[2][0]],
[aij[0][1], aij[0][1] + aij[2][1]],
[aij[0][2], aij[0][2] + aij[2][2]], color = color)
# ay -> ay + ax
ax.plot([aij[1][0], aij[1][0] + aij[0][0]],
[aij[1][1], aij[1][1] + aij[0][1]],
[aij[1][2], aij[1][2] + aij[0][2]], color = color)
# ay -> ay + az
ax.plot([aij[1][0], aij[1][0] + aij[2][0]],
[aij[1][1], aij[1][1] + aij[2][1]],
[aij[1][2], aij[1][2] + aij[2][2]], color = color)
# az -> az + ax
ax.plot([aij[2][0], aij[2][0] + aij[0][0]],
[aij[2][1], aij[2][1] + aij[0][1]],
[aij[2][2], aij[2][2] + aij[0][2]], color = color)
# az -> ax + ay
ax.plot([aij[2][0], aij[2][0] + aij[1][0]],
[aij[2][1], aij[2][1] + aij[1][1]],
[aij[2][2], aij[2][2] + aij[1][2]], color = color)
# ax + ay -> ax + ay + az
ax.plot([aij[0][0] + aij[1][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[0][1] + aij[1][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[0][2] + aij[1][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
# ax + az -> ax + ay + az
ax.plot([aij[0][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[0][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[0][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
# ay + az -> ax + ay + az
ax.plot([aij[1][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[1][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[1][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
[ドキュメント]
def draw_unitcell(ax, sites, aij, nrange, color = 'black'):
"""
概要: 結晶の単位胞とサイト（原子）を3Dプロットに描画します。
詳細説明:
単位胞の枠を描画し、さらに指定された`nrange`の範囲内に存在する
サイト（原子）を球としてプロットします。
サイトの位置は、分数座標から直交座標に変換されます。
:param ax: matplotlib.axes.Axesオブジェクト。3Dプロットを描画する軸。
:type ax: mpl_toolkits.mplot3d.axes3d.Axes3D
:param sites: サイト情報を含むリスト。各サイトは[原子名, サイトラベル, 原子番号, 原子質量, 電荷, 半径, 色, 座標]のリストまたはタプル。Noneの場合、サイトは描画されません。
:type sites: list[list] or None
:param aij: 3x3の行列で、各行が結晶格子ベクトル (ax, ay, az) を表します。
:type aij: numpy.ndarray
:param nrange: 描画する単位胞の範囲を定義するリスト (例: [[min_x, max_x], [min_y, max_y], [min_z, max_z]])。
:type nrange: list
:param color: 単位胞の枠の辺の色。デフォルトは'black'です。
:type color: str, optional
:returns: なし
:rtype: None
"""
draw_box(ax, aij, nrange, color)
if sites is None:
return
for site in sites:
name, label, z, M, q, r, color, pos = 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 posn[0] < nrange[0][0] or nrange[0][1] < posn[0] \
or posn[1] < nrange[1][0] or nrange[1][1] < posn[1] \
or posn[2] < nrange[2][0] or nrange[2][1] < posn[2]:
continue
x, y, z = fractional_to_cartesian(posn, aij)
ax.scatter([x], [y], [z], marker = 'o', c = color, s = kr *r)
[ドキュメント]
def main():
"""
概要: 結晶の単位胞と逆格子空間の単位胞を描画するメイン処理です。
詳細説明:
グローバル変数`lattice_parameters`と`sites`に基づき、
結晶の格子ベクトル、メトリックテンソル、体積を計算し表示します。
その後、実空間の単位胞と原子、およびスケール調整された逆空間の単位胞を
matplotlibの3Dプロットに描画し、結果を表示します。
:returns: なし
:rtype: None
"""
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))
fig = plt.figure(figsize = figsize)
ax = fig.add_subplot(111, projection='3d')
# Real space unit cell
draw_unitcell(ax, sites, aij, nrange)
# 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, color = '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()