interatomic_distance.py ドキュメント
概要
原子間距離を計算するためのスクリプト。
詳細
結晶構造の格子パラメータとサイト情報に基づいて、指定された距離範囲内の原子間距離を計算し、表示します。
tkcrystalbase モジュールに依存しており、結晶学的な計算（格子ベクトル、計量テンソルなど）はそのモジュールが提供する関数を使用します。
計算された原子間距離は、サイトの組み合わせと周期境界条件を考慮した単位胞の並進ベクトルと共にソートして出力されます。
非標準ライブラリ
このプログラムは以下の非標準ライブラリに依存しています。
numpy: 数値計算に使用されます。特に、配列操作、数学関数（sin,cos,tan,arcsin,arccos,arctan,exp,log,sqrt）、線形代数（linalg）機能を利用します。matplotlib: 3Dプロット（mpl_toolkits.mplot3d.Axes3D）を含め、データ可視化に使用されますが、本スクリプトではプロット処理は実装されていません。tkcrystalbase: 結晶学的な計算（例:cal_lattice_vectors,cal_metrics,cal_volume,cal_reciprocal_lattice_vectors,cal_reciprocal_lattice_parameters,distance）を提供するカスタムモジュールです。
グローバル変数および定数
プログラム内で定義されているグローバル変数および定数は以下の通りです。
lattice_parameters: 結晶の格子パラメータを格納するリストです。形式:
[a, b, c, alpha, beta, gamma]単位: オングストローム (
a,b,c) および度 (alpha,beta,gamma)例:
[5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
sites: 結晶構造内の原子サイト情報を格納するリストです。各サイトは以下の情報を含むリストで表現されます:
[原子名 (str), サイトラベル (str), 原子番号 (int), 原子質量 (float), 電荷 (float), 半径 (float), 色 (str), 単位胞内位置 (numpy.array)]
rmin: 原子間距離計算の下限値（オングストローム）。同一サイトの判定にも使用されます。例:
0.1
rmax: 原子間距離計算の上限値（オングストローム）。例:
4.5
関数
main()
概要
プログラムの主要な処理を実行し、原子間距離を計算して出力します。
詳細
この関数は以下のステップを実行します。
与えられた格子パラメータから格子ベクトル、計量テンソル、単位胞体積、およびそれらの逆格子版を計算し、標準出力に表示します。
原子間距離の計算範囲
rmaxに基づいて、考慮すべき単位胞の最大並進数 (nxmax,nymax,nzmax) を決定します。sitesリストに定義された全ての原子サイトのペアについて、計算された単位胞の並進範囲内で原子間距離を計算します。計算された距離が
rminとrmaxの間にある場合、その距離と関連するサイト情報、並進ベクトルをrlistに追加します。rlistを距離でソートした後、結果を整形して標準出力に表示します。
入力
この関数は引数を取りませんが、以下のグローバル変数に定義された値を使用します。
lattice_parameterssitesrminrmax
出力
この関数は、計算結果を標準出力にテキスト形式で出力します。出力される情報の種類は以下の通りです。
格子パラメータと関連情報:
初期設定の格子パラメータ
計算された格子ベクトル (
aij)計量テンソル (
gij)単位胞体積
逆格子パラメータと関連情報:
逆格子パラメータ (
Rlatt)逆格子ベクトル (
Raij)逆格子計量テンソル (
Rgij)逆格子単位胞体積
単位胞並進範囲:
rmaxに基づいて決定された最大並進数 (nxmax,nymax,nzmax)
原子間距離:
rminとrmaxの範囲内にある全ての原子間距離各距離は、以下の形式で表示されます:
サイトラベル0 (x0, y0, z0) - サイトラベル1 (x1, y1, z1) + (並進x, 並進y, 並進z): dis = 距離 A
処理詳細
初期情報の表示と格子計算:
グローバル変数
lattice_parametersを表示します。tkcrystalbase.cal_lattice_vectors()を使用して格子ベクトルaijを計算し、標準出力に表示します。tkcrystalbase.cal_metrics()を使用して計量テンソルgijを計算し、標準出力に表示します。tkcrystalbase.cal_volume()を使用して単位胞体積を計算し、標準出力に表示します。
逆格子情報の計算と表示:
tkcrystalbase.cal_reciprocal_lattice_vectors()を使用して逆格子ベクトルRaijを計算し、標準出力に表示します。tkcrystalbase.cal_reciprocal_lattice_parameters()を使用して逆格子パラメータRlattを計算し、標準出力に表示します。tkcrystalbase.cal_metrics()を使用して逆格子計量テンソルRgijを計算し、標準出力に表示します。tkcrystalbase.cal_volume()を使用して逆格子単位胞体積を計算し、標準出力に表示します。
単位胞並進範囲の決定:
rmaxとlattice_parametersの値に基づいて、X, Y, Z方向の最大並進数 (nxmax,nymax,nzmax) を計算し、標準出力に表示します。
原子間距離の計算:
sitesリスト内の全てのサイトペア (site0とsite1) についてループします。各サイトペアに対し、-
nmaxから+nmaxまでの並進ベクトル (ix,iy,iz) を考慮してループします。tkcrystalbase.distance()関数を使用して、並進されたサイト間の距離を計算します。計算された距離が
rminより大きくrmaxより小さい場合、その距離、両サイトの情報、および並進ベクトルをrlistリストに追加します。
結果のソートと出力:
rlistを距離 (x[0]) を主キーとし、次にサイト情報、並進ベクトルでソートします。ソートされた
rlistの各エントリについて、整形された原子間距離情報を標準出力に表示します。プログラムは
exit()を呼び出して終了します。
使用例
このスクリプトは、Pythonインタープリタで直接実行することを想定しています。
python interatomic_distance.py
実行すると、計算された格子情報、逆格子情報、および原子間距離が標準出力に表示されます。
ソースコード
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 *
"""
概要: 原子間距離を計算するためのスクリプト。
詳細説明:
結晶構造の格子パラメータとサイト情報に基づいて、指定された距離範囲内の原子間距離を計算し、表示します。
TKCrystalBaseモジュールに依存しており、結晶学的な計算（格子ベクトル、計量テンソルなど）はそのモジュールが提供する関数を使用します。
計算された原子間距離は、サイトの組み合わせと周期境界条件を考慮した単位胞の並進ベクトルと共にソートして出力されます。
関連リンク: :doc:`interatomic_distance_usage`
"""
# 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])]
]
# Distance range
rmin = 0.1 # angstrom. used to judge identical site
rmax = 4.5 # angstrom.
def main():
"""
概要: プログラムの主要な処理を実行し、原子間距離を計算して出力する。
詳細説明:
この関数は以下のステップを実行します。
1. 与えられた格子パラメータから格子ベクトル、計量テンソル、単位胞体積、およびそれらの逆格子版を計算し、標準出力に表示します。
2. 原子間距離の計算範囲 `rmax` に基づいて、考慮すべき単位胞の最大並進数 (`nxmax`, `nymax`, `nzmax`) を決定します。
3. `sites` リストに定義された全ての原子サイトのペアについて、計算された単位胞の並進範囲内で原子間距離を計算します。
4. 計算された距離が `rmin` と `rmax` の間にある場合、その距離と関連するサイト情報、並進ベクトルを `rlist` に追加します。
5. `rlist` を距離でソートした後、結果を整形して標準出力に表示します。
"""
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))
# Calculate the range of unit cells
nxmax = int(rmax / lattice_parameters[0]) + 1
nymax = int(rmax / lattice_parameters[1]) + 1
nzmax = int(rmax / lattice_parameters[2]) + 1
print("")
print("nmax:", nxmax, nymax, nzmax)
# Calculate interatomic distances and store them in rlist list variable
rlist = []
for isite0 in range(len(sites)):
site0 = sites[isite0]
name0, label0, z0, M0, q0, r0, color0, pos0 = site0
for isite1 in range(len(sites)):
site1 = sites[isite1]
name1, label1, z1, M1, q1, r1, color1, pos1 = site1
for iz in range(-nzmax, nzmax+1):
for iy in range(-nymax, nymax+1):
for ix in range(-nxmax, nxmax+1):
dis = distance(pos0, pos1 + np.array([ix, iy, iz]), gij)
if dis < rmin or rmax < dis:
continue
rlist.append([dis, site0, site1, ix, iy, iz])
# print(" {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
# .format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))
# Sort rlist by distance (x[0] priority)
rlist.sort(key = lambda x: (x[0], x[1], x[3], x[4], x[5], x[2]))
print("")
print("Interatomic distances:")
for rinf in rlist:
dis, site0, site1, ix, iy, iz = rinf
name0, label0, z0, M0, q0, r0, color0, pos0 = site0
name1, label1, z1, M1, q1, r1, color1, pos1 = site1
print(" {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
.format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))
print("")
exit()
if __name__ == '__main__':
main()