"""
結晶学計算のための基本ユーティリティ関数と定数を提供するモジュール。
本モジュールには、物理定数、格子パラメータ、サイト情報、
および格子ベクトル、逆格子ベクトル、距離、角度などの計算関数が含まれます。
: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未満の範囲に削減する。
引数からその整数部分を引くことで、実数を[0, 1)の範囲にマッピングします。
:param x: 削減する数値。
:type x: float
:returns: 0以上1未満に削減された数値。
:rtype: float
"""
return x - int(x)
return x
[ドキュメント]
def round01(x):
"""
数値を0または1に丸める、または0以上1未満の範囲に削減する。
許容誤差に基づいて0または1に丸めます。
1.0に近い場合は1.0に、0.0に近い場合は0.0に丸めます。
それ以外の場合、1.0より大きい場合は整数部分を引いて削減し、
1.0より小さい場合は整数部分を引いた後に1.0を足して0以上1未満の範囲に調整します。
:param x: 丸める数値。
:type x: float
:returns: 丸められた、または削減された数値。
:rtype: float
"""
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):
"""
数値を指定された許容誤差の倍数に丸める。
指定された許容誤差 `tol` の倍数に `x` を丸めます。
丸めの際には `0.1 * tol` が加算されるため、丸めの基準がわずかにシフトします。
:param x: 丸める数値。
:type x: float
:param tol: 許容誤差（丸めの単位）。
:type tol: float
:returns: 許容誤差の倍数に丸められた数値。
:rtype: float
"""
val = tol * int( (x+0.1*tol) / tol )
return val
[ドキュメント]
def cal_lattice_vectors(latt):
"""
格子定数から実空間の格子ベクトルを計算する。
六方晶系または三斜晶系の格子定数 (a, b, c, alpha, beta, gamma) から、
結晶学的な格子ベクトル `aij` (a1, a2, a3) を計算します。
格子定数はアンゴストローム、角度は度単位で与えられます。
:param latt: [a, b, c, alpha, beta, gamma] の形式の格子定数リスト。角度は度単位。
:type latt: list[float]
:returns: 3x3のnumpy配列として表された格子ベクトル。各行がベクトルa1, a2, a3に対応します。
:rtype: numpy.ndarray
"""
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] * cosb - aij[2][0] * cosg)
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):
"""
格子定数から格子ベクトルと計量テンソルを計算する。
`cal_lattice_vectors` を使用して格子ベクトルを計算し、
それらから計量テンソル `gij` を導出します。
結果は辞書形式で返されます。
:param latt: [a, b, c, alpha, beta, gamma] の形式の格子定数リスト。角度は度単位。
:type latt: list[float]
:returns: 格子ベクトル (aij) と計量テンソル (gij) を含む辞書。
:rtype: dict
* `'aij'` (`numpy.ndarray`): 格子ベクトル (3x3配列)。
* `'gij'` (`numpy.ndarray`): 計量テンソル (3x3配列)。
"""
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つの格子ベクトル (a1, a2, a3) からなる行列の混合積 (a1 . (a2 x a3)) を計算して体積を求めます。
:param aij: 3x3のnumpy配列として表された格子ベクトル。各行がベクトルa1, a2, a3に対応します。
:type aij: numpy.ndarray
:returns: 単位胞の体積。
:rtype: 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):
"""
実空間の格子ベクトルから逆格子ベクトルを計算する。
以下の式に基づき逆格子ベクトル `Ra`, `Rb`, `Rc` を計算します。
Ra = (a2 x a3) / V
Rb = (a3 x a1) / V
Rc = (a1 x a2) / V
ここで `V` は実空間の単位胞体積です。
:param aij: 3x3のnumpy配列として表された実空間の格子ベクトル。
:type aij: numpy.ndarray
:returns: 逆格子ベクトル [Ra, Rb, Rc] のリスト。各要素はnumpy配列。
:rtype: list[numpy.ndarray]
"""
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):
"""
逆格子ベクトルから逆格子定数を計算する。
逆格子ベクトル `Ra`, `Rb`, `Rc` から、
逆格子定数 (`Ra`, `Rb`, `Rc`, `Ralpha`, `Rbeta`, `Rgamma`) を計算します。
角度は度単位で返されます。
:param Raij: 逆格子ベクトル [Ra, Rb, Rc] のリスト。
:type Raij: list[numpy.ndarray]
:returns: [Ra, Rb, Rc, Ralpha, Rbeta, Rgamma] の形式の逆格子定数リスト。角度は度単位。
:rtype: list[float]
"""
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):
"""
分数座標をデカルト座標に変換する。
単位胞内の分数座標 `(u, v, w)` を、与えられた格子ベクトル `aij` を用いて
デカルト座標 `(x, y, z)` に変換します。
x = u * a1x + v * a2x + w * a3x
y = u * a1y + v * a2y + w * a3y
z = u * a1z + v * a2z + w * a3z
:param pos: 変換する分数座標 (u, v, w)。
:type pos: numpy.ndarray
:param aij: 3x3のnumpy配列として表された格子ベクトル。
:type aij: numpy.ndarray
:returns: デカルト座標 (x, y, z)。
:rtype: tuple[float, float, float]
"""
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つの分数座標間の距離の二乗を計算する。
計量テンソル `gij` を用いて、2つの分数座標 `pos0` と `pos1` 間の
ユークリッド距離の二乗を計算します。
:param pos0: 1点目の分数座標。
:type pos0: numpy.ndarray
:param pos1: 2点目の分数座標。
:type pos1: numpy.ndarray
:param gij: 3x3のnumpy配列として表された計量テンソル。
:type gij: numpy.ndarray
:returns: 2点間の距離の二乗。
:rtype: float
"""
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つの分数座標間の距離を計算する。
`distance2` 関数で計算された距離の二乗の平方根を取ることで、2点間の距離を計算します。
:param pos0: 1点目の分数座標。
:type pos0: numpy.ndarray
:param pos1: 2点目の分数座標。
:type pos1: numpy.ndarray
:param gij: 3x3のnumpy配列として表された計量テンソル。
:type gij: numpy.ndarray
:returns: 2点間の距離。
:rtype: float
"""
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` を頂点として `pos1` と `pos2` がなす角度を、
計量テンソル `gij` を用いて計算します。
頂点からいずれかの点までの距離が0の場合は0.0を返します。
計算された角度は度単位で返されます。
:param pos0: 角度の頂点となる分数座標。
:type pos0: numpy.ndarray
:param pos1: 頂点から伸びる1つ目のベクトルの終点の分数座標。
:type pos1: numpy.ndarray
:param pos2: 頂点から伸びる2つ目のベクトルの終点の分数座標。
:type pos2: numpy.ndarray
:param gij: 3x3のnumpy配列として表された計量テンソル。
:type gij: numpy.ndarray
:returns: 3点がなす角度（度単位）。
:rtype: float
"""
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;
[ドキュメント]
def main():
"""
モジュールのサンプル実行関数。
定義されている格子パラメータとサイト情報を用いて、
格子ベクトル、計量テンソル、体積、逆格子関連の計算結果、
および距離と角度の計算例を出力します。
最後にプログラムを終了します。
:returns: None
: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]))
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()