"""
tkcrystalbase.py
概要: 結晶構造計算と描画に関連するユーティリティ関数を提供するモジュール。
詳細説明:
このモジュールは、結晶の格子定数から実空間および逆空間の格子ベクトル、
メトリックテンソル、単位格子の体積などを計算する機能を提供します。
また、原子の分数座標とデカルト座標間の変換、原子間距離や角度の計算、
そしてmatplotlibを使用した3D描画ヘルパー関数を含んでいます。
物理定数や一般的な格子パラメータ、サイト情報もここで定義されています。
関連リンク:
:doc:`tkcrystalbase_usage`
"""
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from pprint import pprint
pi = 3.14159265358979323846
pi2 = pi + pi
torad = 0.01745329251944 # rad/deg";
todeg = 57.29577951472 # deg/rad";
basee = 2.71828183
h = 6.6260755e-34 # Js";
h_bar = 1.05459e-34 # "Js";
hbar = h_bar
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
me = 9.1093897e-31 # kg";
mp = 1.6726231e-27 # kg";
mn = 1.67495e-27 # kg";
u0 = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
e2_4pie0 = 2.30711e-28 # Nm<sup>2</sup>";
a0 = 5.29177e-11 # m";
kB = 1.380658e-23 # JK<sup>-1</sup>";
NA = 6.0221367e23 # mol<sup>-1</sup>";
R = 8.31451 # J/K/mol";
F = 96485.3 # C/mol";
g = 9.81 # m/s2";
# 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])]
]
[ドキュメント]
def reduce01(x):
"""
概要: 数値を0以上1未満の範囲に変換する。
詳細説明:
入力された数値の小数部分を返します。
例えば、3.75は0.75に、-2.3は0.7に変換されます。
:param x: float: 変換する数値。
:returns: float: 0以上1未満の範囲に変換された値。
"""
return x - int(x)
return x
[ドキュメント]
def round01(x):
"""
概要: 数値を特定の閾値に基づいて0または1に丸めるか、0以上1未満の範囲に変換する。
詳細説明:
`x` が1.0に近い（絶対差が0.0002未満）場合は1.0を返します。
`x` が0.0に近い（絶対差が0.0002未満）場合は0.0を返します。
それ以外の場合、`x > 1.0` なら `x - int(x)` を返します。
`x < 1.0` なら `x - int(x) + 1.0` を返します。
この関数は、負の数や1を超える数値を0以上1未満の範囲に正規化するために使用されることがあります。
:param x: float: 丸める、または変換する数値。
:returns: float: 丸められた、または0以上1未満の範囲に変換された値。
"""
if abs(x - 1.0) < 0.0002:
return 1.0
if abs(x) < 0.0002:
return 0.0
if x > 1.0:
return x - int(x)
if x < 1.0:
return x - int(x) + 1.0
return x
[ドキュメント]
def round_parameter(x, tol):
"""
概要: 数値を指定された許容範囲 (tolerance) の整数倍に丸める。
詳細説明:
入力値 `x` を `tol` の整数倍に丸めます。
丸め処理は、`(x + 0.1 * tol) / tol` を整数に変換し、`tol` を掛けることで行われます。
これにより、例えば `tol=0.5` で `x=1.2` なら `1.0` に、`x=1.3` なら `1.5` に丸められます。
:param x: float: 丸める数値。
:param tol: float: 丸めの基準となる許容範囲（単位）。
:returns: float: `tol` の整数倍に丸められた値。
"""
val = tol * int( (x+0.1*tol) / tol )
return val
[ドキュメント]
def cal_lattice_vectors(latt):
"""
概要: 単位格子の格子定数から実空間の格子ベクトルを計算する。
詳細説明:
格子定数 (a, b, c, α, β, γ) を受け取り、
直交座標系における3つの格子ベクトル a_vec, b_vec, c_vec を含む3x3行列を計算します。
格子ベクトルは、x軸に沿って a_vec が配置されるように構成されます。
角度は度数で指定されますが、内部ではラジアンに変換して計算されます。
:param latt: list of float: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
長さはオングストローム、角度は度で与えられます。
:returns: numpy.ndarray: 3x3 の格子ベクトル行列。
`[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
"""
cosa = cos(torad * latt[3])
cosb = cos(torad * latt[4])
cosg = cos(torad * latt[5])
sing = sin(torad * latt[5])
aij = np.empty([3, 3], dtype = float)
aij[0][0] = latt[0]
aij[0][1] = 0.0;
aij[0][2] = 0.0;
aij[1][0] = latt[1] * cosg
aij[1][1] = latt[1] * sing
aij[1][2] = 0.0;
aij[2][0] = latt[2] * cosb
aij[2][1] = latt[2] * (cosa - cosb * cosg) / sing
if abs(aij[2][1]) < 1.0e-8:
aij[2][1] = 0.0
else:
aij[2][1] = aij[2][1] / sing
aij[2][2] = sqrt(latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1])
return aij
[ドキュメント]
def cal_metrics(latt):
"""
概要: 格子定数から格子ベクトルとメトリックテンソルを計算する。
詳細説明:
単位格子の格子定数から実空間の格子ベクトル `aij` を計算し、
それらの内積からメトリックテンソル `gij` を構築します。
結果は辞書形式で返されます。
:param latt: list of float: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
長さはオングストローム、角度は度で与えられます。
:returns: dict: 計算された格子ベクトル行列 (`'aij'`) とメトリックテンソル (`'gij'`)
をキーに持つ辞書。
`{'aij': numpy.ndarray, 'gij': numpy.ndarray}` の形式。
"""
inf = {}
cosa = cos(torad * latt[3])
cosb = cos(torad * latt[4])
cosg = cos(torad * latt[5])
aij = cal_lattice_vectors(latt)
inf['aij'] = aij
gij = np.empty([3, 3], dtype = float)
for i in range(3):
for j in range(i, 3):
gij[i][j] = np.dot(aij[i], aij[j])
gij[j][i] = gij[i][j]
inf['gij'] = gij
return inf
[ドキュメント]
def cal_volume(aij):
"""
概要: 格子ベクトルから単位格子の体積を計算する。
詳細説明:
3つの格子ベクトル (a_vec, b_vec, c_vec) のスカラー三重積 (a_vec・(b_vec × c_vec)) を用いて、
単位格子の体積を計算します。
:param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
`[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
:returns: float: 単位格子の体積。
"""
axb = np.cross(aij[0], aij[1]) # Outner product
vol = np.dot(axb, aij[2]) # Inner product
return vol
[ドキュメント]
def cal_reciprocal_lattice_vectors(aij):
"""
概要: 実空間格子ベクトルから逆空間格子ベクトルを計算する。
詳細説明:
実空間の格子ベクトル `aij` (a_vec, b_vec, c_vec) を用いて、
逆空間の格子ベクトル (a*_vec, b*_vec, c*_vec) を計算します。
逆空間ベクトルは、単位格子の体積 V と実空間ベクトルの外積によって定義されます。
:param aij: numpy.ndarray: 3x3 の実空間格子ベクトル行列。
`[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
:returns: list of numpy.ndarray: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
"""
V = cal_volume(aij)
Ra = np.cross(aij[1], aij[2]) / V
Rb = np.cross(aij[2], aij[0]) / V
Rc = np.cross(aij[0], aij[1]) / V
return [Ra, Rb, Rc]
[ドキュメント]
def cal_reciprocal_lattice_parameters(Raij):
"""
概要: 逆空間格子ベクトルから逆空間の格子定数を計算する。
詳細説明:
逆空間の格子ベクトル `Raij` (a*_vec, b*_vec, c*_vec) を用いて、
逆空間の格子定数 (a*, b*, c*, α*, β*, γ*) を計算します。
長さはベクトルのノルムから、角度はベクトルの内積から導出されます。
:param Raij: list of numpy.ndarray: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
:returns: list of float: 逆空間の格子定数 `[Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]` のリスト。
"""
Ra = la.norm(Raij[0])
Rb = la.norm(Raij[1])
Rc = la.norm(Raij[2])
Ralpha = todeg * arccos(np.dot(Raij[1], Raij[2]) / Rb / Rc)
Rbeta = todeg * arccos(np.dot(Raij[2], Raij[0]) / Rc / Ra)
Rgamma = todeg * arccos(np.dot(Raij[0], Raij[1]) / Ra / Rb)
return [Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]
[ドキュメント]
def fractional_to_cartesian(pos, aij):
"""
概要: 分数座標からデカルト座標に変換する。
詳細説明:
単位格子内の原子の分数座標 (x_frac, y_frac, z_frac) と格子ベクトル `aij` を用いて、
原子のデカルト座標 (x_cart, y_cart, z_cart) を計算します。
:param pos: numpy.ndarray or list of float: 変換する分数座標 `[x_frac, y_frac, z_frac]`。
:param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
`[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
:returns: tuple of float: 計算されたデカルト座標 `(x_cart, y_cart, z_cart)`。
"""
x = pos[0] * aij[0][0] + pos[1] * aij[1][0] + pos[2] * aij[2][0]
y = pos[0] * aij[0][1] + pos[1] * aij[1][1] + pos[2] * aij[2][1]
z = pos[0] * aij[0][2] + pos[1] * aij[1][2] + pos[2] * aij[2][2]
return x, y, z
[ドキュメント]
def distance2(pos0, pos1, gij):
"""
概要: 2つの分数座標間の距離の2乗を計算する。
詳細説明:
2つの原子の分数座標 `pos0` と `pos1` の間の距離の2乗を、
メトリックテンソル `gij` を使用して計算します。
この関数は周期境界条件を考慮しません。
:param pos0: numpy.ndarray or list of float: 1点目の分数座標 `[x0, y0, z0]`。
:param pos1: numpy.ndarray or list of float: 2点目の分数座標 `[x1, y1, z1]`。
:param gij: numpy.ndarray: 3x3 のメトリックテンソル。
:returns: float: 2点間の距離の2乗。
"""
dx = pos1 - pos0
# dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
+ 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
return r2
[ドキュメント]
def distance(pos0, pos1, gij):
"""
概要: 2つの分数座標間の距離を計算する。
詳細説明:
2つの原子の分数座標 `pos0` と `pos1` の間の距離を、
メトリックテンソル `gij` を使用して計算します。
この関数は周期境界条件を考慮しません。
:param pos0: numpy.ndarray or list of float: 1点目の分数座標 `[x0, y0, z0]`。
:param pos1: numpy.ndarray or list of float: 2点目の分数座標 `[x1, y1, z1]`。
:param gij: numpy.ndarray: 3x3 のメトリックテンソル。
:returns: float: 2点間の距離。
"""
dx = pos1 - pos0
# dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
+ 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
r = sqrt(r2)
return r
[ドキュメント]
def angle(pos0, pos1, pos2, gij):
"""
概要: 3つの分数座標で定義される角度を計算する。
詳細説明:
`pos0` を頂点とし、`pos0-pos1` と `pos0-pos2` の2つのベクトルがなす角度を計算します。
角度の計算にはメトリックテンソル `gij` が使用されます。
結果は度数で返され、180度を超える場合は `360.0 - angle` として調整されます。
:param pos0: numpy.ndarray or list of float: 角度の頂点となる分数座標 `[x0, y0, z0]`。
:param pos1: numpy.ndarray or list of float: 1つ目のベクトルを定義する分数座標 `[x1, y1, z1]`。
:param pos2: numpy.ndarray or list of float: 2つ目のベクトルを定義する分数座標 `[x2, y2, z2]`。
:param gij: numpy.ndarray: 3x3 のメトリックテンソル。
:returns: float: 3点間で形成される角度（度）。
いずれかの距離が0の場合は0.0を返します。
"""
dis01 = distance(pos0, pos1, gij)
if dis01 == 0.0:
return 0.0
dis02 = distance(pos0, pos2, gij)
if dis02 == 0.0:
return 0.0
dx01 = pos1 - pos0
dx02 = pos2 - pos0
# dx01 = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
# dx02 = [pos2[0] - pos0[0], pos2[1] - pos0[1], pos2[2] - pos0[2]]
ip = gij[0][0] * dx01[0]*dx02[0] + gij[1][1] * dx01[1]*dx02[1] + gij[2][2] * dx01[2]*dx02[2] \
+ 2.0 * (gij[0][1] * dx01[0]*dx02[1] + gij[0][2] * dx01[0]*dx02[2] + gij[1][2] * dx01[1]*dx02[2])
cosa = ip / dis01 / dis02
angle = todeg * arccos(cosa)
if angle > 180.0:
angle = 360.0 - angle
return angle;
# ax.set_box_aspect([xrange[1] - xrange[0], yrange[1] - yrange[0], zrange[1] - zrange[0]])
[ドキュメント]
def draw_box(ax, aij, nrange, color = 'black'):
"""
概要: 3D軸上に単位格子のボックスを描画する。
詳細説明:
格子ベクトル `aij` で定義される単位格子の12本の辺を、Matplotlibの3D軸 `ax` に描画します。
この関数は、単一の単位格子を描画するために設計されており、`nrange` は現在の実装では使用されません。
:param ax: matplotlib.axes.Axes: 描画対象のMatplotlib 3D軸オブジェクト。
:param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
`[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
:param nrange: list of list of int: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
この関数では使用されません。
:param color: str, optional: ボックスの線の色。デフォルトは 'black'。
:returns: 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', facecolor = 'black', edgecolor = 'white', alpha = 0.7, kr = 1.0):
"""
概要: 3D軸上に単位格子とその内部のサイト原子を描画する。
詳細説明:
まず `draw_box` 関数を呼び出して単位格子の境界線を描画します。
次に、`sites` リストに含まれる各原子について、その分数座標をデカルト座標に変換し、
指定された `nrange` の範囲内の周期的な位置に原子を散布図として描画します。
原子の描画スタイル（色、透明度、サイズ）は引数でカスタマイズ可能です。
:param ax: matplotlib.axes.Axes: 描画対象のMatplotlib 3D軸オブジェクト。
:param sites: list of list: サイト情報（原子名、ラベル、原子番号、原子質量、電荷、半径、色、位置）のリスト。
各要素は `[name, label, z, M, q, r, color, pos]` の形式。
`None` の場合、原子は描画されません。
:param aij: numpy.ndarray: 3x3 の格子ベクトル行列。
`[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
:param nrange: list of list of int: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
原子の繰り返し構造を描画するために使用されます。
:param color: str, optional: 単位格子ボックスの線の色。デフォルトは 'black'。
:param facecolor: str, optional: サイト原子の面の描画色。デフォルトは 'black'。
:param edgecolor: str, optional: サイト原子のエッジの描画色。デフォルトは 'white'。
:param alpha: float, optional: サイト原子の透明度（0.0から1.0）。デフォルトは 0.7。
:param kr: float, optional: サイト原子の描画サイズに対するスケール因子。デフォルトは 1.0。
:returns: None
"""
draw_box(ax, aij, nrange, color)
if sites is None: return
for isite, site in enumerate(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 = facecolor, edgecolors = edgecolor, alpha = alpha, s = kr *r)
[ドキュメント]
def main():
"""
概要: モジュールの主要な計算と結果の表示を実行する。
詳細説明:
定義された `lattice_parameters` を使用して、実空間および逆空間の
格子ベクトル、メトリックテンソル、単位格子の体積を計算します。
さらに、サンプルとして原子間距離と角度の計算も行い、
その結果を標準出力に表示します。
最後にプログラムを終了します。
: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]))
print("")
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))
print("")
print("dis=", distance(np.array([0,0,0]), np.array([1,1,1]), gij))
print("angle=", angle (np.array([0,0,0]), np.array([1,1,1]), np.array([1,0,0]), gij))
print("angle=", angle (np.array([0,0,0]), np.array([1,0,0]), np.array([0,1,0]), gij))
print("")
exit()
if __name__ == '__main__':
main()