crystal_distance.py ダウンロード/コピー
crystal_distance.py
crystal_distance.py
1import sys
2import os
3from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
4import numpy as np
5from numpy import linalg as la
6from mpl_toolkits.mplot3d import Axes3D
7import matplotlib.pyplot as plt
8
9from tkcrystalbase import *
10
11
12"""
13結晶中の原子間距離を計算するスクリプト。
14
15このスクリプトは、与えられた格子定数と原子サイト情報に基づき、結晶中の原子間距離を計算し、表示します。
16実空間および逆空間の格子情報（格子ベクトル、メトリックテンソル、単位胞体積）も計算し、表示します。
17計算された原子間距離は、指定された距離範囲内でソートされて出力されます。
18
19:doc:`crystal_distance_usage`
20"""
21
22
23# Lattice parameters (angstrom and degree)
24#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
25lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
26
27# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
28sites = [
29 ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
30 ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.5, 0.5])]
31 ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.0, 0.5])]
32 ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.5, 0.0])]
33 ,['Cl', 'Cl1', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
34 ,['Cl', 'Cl2', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
35 ,['Cl', 'Cl3', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
36 ,['Cl', 'Cl4', 17, 35.4527, -1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
37 ]
38
39# Distance range
40rmin = 0.1 # angstrom. used to judge identical site
41rmax = 4.5 # angstrom.
42
43
44def main():
45 """
46 結晶の格子情報と原子間距離を計算し、結果を出力するメイン関数。
47
48 この関数は、`lattice_parameters` と `sites` グローバル変数で定義された情報に基づき、
49 以下の処理を実行します。
50
51 1. 実空間の格子ベクトル、メトリックテンソル、単位胞体積を計算し、標準出力に表示します。
52 2. 逆空間の格子ベクトル、逆格子メトリックテンソル、逆格子単位胞体積を計算し、標準出力に表示します。
53 3. `rmin` と `rmax` で指定された範囲内で、全ての原子ペア間の原子間距離を計算します。
54 この際、周囲の単位胞（`nxmax`, `nymax`, `nzmax` の範囲）に存在する原子も考慮します。
55 4. 計算された原子間距離を昇順にソートし、各距離と関連するサイト情報、単位胞のオフセットを標準出力に表示します。
56 """
57 print("")
58 print("Lattice parameters:", lattice_parameters)
59 aij = cal_lattice_vectors(lattice_parameters)
60 print("Lattice vectors:")
61 print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
62 print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
63 print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
64 inf = cal_metrics(lattice_parameters)
65 gij = inf['gij']
66 print("Metric tensor:")
67 print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
68 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
69 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
70 print("")
71 volume = cal_volume(aij)
72 print("Volume: {:12.4g} A^3".format(volume))
73
74 print("")
75 print("Unit cell volume: {:12.4g} A^3".format(volume))
76 Raij = cal_reciprocal_lattice_vectors(aij)
77 Rlatt = cal_reciprocal_lattice_parameters(Raij)
78 Rinf = cal_metrics(Rlatt)
79 Rgij = Rinf['gij']
80 print("Reciprocal lattice parameters:", Rlatt)
81 print("Reciprocal lattice vectors:")
82 print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
83 print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
84 print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
85 print("Reciprocal lattice metric tensor:")
86 print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
87 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
88 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
89 Rvolume = cal_volume(Raij)
90 print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
91
92# Calculate the range of unit cells
93 nxmax = int(rmax / lattice_parameters[0]) + 1
94 nymax = int(rmax / lattice_parameters[1]) + 1
95 nzmax = int(rmax / lattice_parameters[2]) + 1
96 print("")
97 print("nmax:", nxmax, nymax, nzmax)
98
99# Calculate interatomic distances and store them in rlist list variable
100 rlist = []
101 for isite0 in range(len(sites)):
102 site0 = sites[isite0]
103 name0, label0, z0, M0, q0, r0, color0, pos0 = site0
104 for isite1 in range(len(sites)):
105 site1 = sites[isite1]
106 name1, label1, z1, M1, q1, r1, color1, pos1 = site1
107 for iz in range(-nzmax, nzmax+1):
108 for iy in range(-nymax, nymax+1):
109 for ix in range(-nxmax, nxmax+1):
110 dis = distance(pos0, pos1 + np.array([ix, iy, iz]), gij)
111 if dis < rmin or rmax < dis:
112 continue
113
114 rlist.append([dis, site0, site1, ix, iy, iz])
115
116# print(" {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
117# .format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))
118
119# Sort rlist by distance (x[0] priority)
120 rlist.sort(key = lambda x: (x[0], x[1], x[3], x[4], x[5], x[2]))
121
122 print("")
123 print("Interatomic distances:")
124 for rinf in rlist:
125 dis, site0, site1, ix, iy, iz = rinf
126 name0, label0, z0, M0, q0, r0, color0, pos0 = site0
127 name1, label1, z1, M1, q1, r1, color1, pos1 = site1
128 print(" {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
129 .format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))
130
131 print("")
132 exit()
133
134
135if __name__ == '__main__':
136 main()