"""
X線回折の計算と結晶情報の表示を行うスクリプト。
概要:
結晶格子パラメータから格子ベクトル、逆格子ベクトル、体積などを計算し、
指定されたX線波長と最大2θ角度に基づいて、可能な回折指数(hkl)と
それに対応する回折角(2θ)および格子面間隔(d)を計算し、出力します。
`tkcrystalbase` モジュールから基本的な結晶学ユーティリティをインポートして使用します。
関連リンク:
:doc:`crystal_XRD_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])]
]
# X-ray wavelength
wl = 1.5405 # angstrom
# G min to remove the origin of reciprocal space
Gmin = 1.0e-5
# 2Theta max
Q2max = 150.0 # degree in 2Theta
[ドキュメント]
def main():
"""
プログラムの主要な処理を実行するメイン関数。
概要:
結晶構造情報に基づいてX線回折の計算を行い、結果を標準出力に表示します。
詳細説明:
結晶格子パラメータ、サイト情報、X線波長、最大2θ角などのグローバル変数を参照し、
以下の処理を行います。
1. 格子ベクトル、メートルテンソル、体積を計算し表示します。
2. 逆格子ベクトル、逆格子パラメータ、逆格子メートルテンソル、逆格子体積を計算し表示します。
3. 指定された最大2θ角から最小格子面間隔(dmin)を計算し、それに基づいてhkl指数の最大範囲を決定します。
4. 計算されたhkl範囲内で、各(h, k, l)に対する逆格子空間での距離(G)、格子面間隔(d)、
X線回折角(2θ)を計算します。
5. 有効な回折ピーク（Gminより大きいG値、sinQが1.0未満、2θmaxより小さい2θ）のみをリストに格納します。
6. 回折ピークのリストを2θ値でソートし、結果を表示します。
:param None:
:returns: 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))
dmin = wl / 2.0 / sin(0.5 * Q2max * torad)
hmax = int(lattice_parameters[0] / dmin)
kmax = int(lattice_parameters[1] / dmin)
lmax = int(lattice_parameters[2] / dmin)
print("")
print("hkl range:", hmax, kmax, lmax)
# Calculate diffraction angles and store them in qlist list variable
org = np.array([0.0, 0.0, 0.0])
qlist = []
for l in range(-lmax, lmax+1):
for k in range(-kmax, kmax+1):
for h in range(-hmax, hmax+1):
# Calculate distance in reciprocal space between (0, 0, 0) and (h, k, l)
G = distance(org, np.array([h, k, l]), Rgij)
if G < Gmin:
continue
# Calculate lattice space from G
d = 1.0 / G
sinQ = wl / 2.0 / d
if sinQ >= 1.0:
continue
# Calculate diffraction angle 2Theta
Q2 = 2.0 * todeg * arcsin(sinQ)
if Q2 > Q2max:
continue
qlist.append([Q2, d, h, k, l])
# print(" 2Q={:12.4g} d={:12.6g} ({:3d} {:3d} {:3d})".format(Q2, d, h, k, l))
# Sort rlist by 2Theta (x[0] priority)
qlist.sort(key = lambda x: (x[0], x[2], x[3], x[4]))
print("")
print("Diffraction angle, d, h, k, l:")
for qinf in qlist:
Q2, d, h, k, l = qinf
print(" 2Q={:12.4g} d={:12.6g} ({:3d} {:3d} {:3d})".format(Q2, d, h, k, l))
print("")
exit()
if __name__ == '__main__':
main()