symmetrize_tklib.py ダウンロード/コピー
symmetrize_tklib.py
symmetrize_tklib.py
1"""
2pymatgenとtklibを用いてCIFファイルを対称化し、その結果を新しいCIFファイルに保存するスクリプト。
3
4概要:
5 入力されたCIFファイルを読み込み、pymatgenのSpacegroupAnalyzerを使用して結晶構造を対称化します。
6 対称化された構造から得られる空間群情報、格子定数、原子サイトの情報を抽出し、
7 tklibのtkFileモジュールを使用して新しいCIFファイルとして出力します。
8 このスクリプトは、pymatgenの強力な対称性解析機能を活用しつつ、tkCrytalのデータ構造と連携して
9 出力を生成するハイブリッドなアプローチを採用しています。
10
11詳細説明:
12 1. コマンドライン引数から入力CIFファイルを指定します（デフォルトは 'SrTiO3.cif'）。
13 2. ログファイルを生成し、標準出力とログファイルへ出力をリダイレクトします。
14 3. tkCIFを使って入力CIFファイルを読み込み、基本的な結晶情報を取得します。
15 4. pymatgenのStructure.from_fileメソッドで入力CIFファイルを読み込み、Structureオブジェクトを生成します。
16 5. SpacegroupAnalyzerを用いてStructureオブジェクトを解析し、対称化されたStructureオブジェクトを取得します。
17 6. 対称化された構造から、空間群名、番号、対称操作、格子パラメータ、等価サイトの情報を抽出します。
18 7. 抽出した情報（格子、空間群、対称操作、原子サイトの分数座標と占有率）を基に、
19 新しいCIFファイル（'{filebody}-symmetrized.cif'）を生成・保存します。
20 8. 対称操作のマトリクス成分をCIF形式の文字列（例: 'x', 'x+y', '1/2-z'）に変換するヘルパー関数を使用します。
21
22関連リンク:
23 :doc:`symmetrize_tklib_usage`
24"""
25import re
26import sys
27
28from pymatgen.io.cif import CifParser
29from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
30from pymatgen.core.structure import Structure
31from pymatgen.core.lattice import Lattice
32from pymatgen.core.periodic_table import Element
33
34
35from tklib.tkapplication import tkApplication
36from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
37from tklib.tkfile import tkFile
38from tklib.tkcrystal.tkcif import tkCIF
39
40
41infile = 'SrTiO3.cif'
42prec = 1.0e-3
43
44app = tkApplication()
45argv = sys.argv
46narg = len(argv)
47if narg >= 2:
48 infile = argv[1]
49
50
51#==========================================
52# Main prgram
53#==========================================
54def main():
55 """
56 メインプログラムのエントリーポイント。CIFファイルを読み込み、対称化し、結果を出力する。
57
58 詳細説明:
59 - 入力CIFファイルを読み込み、`tkCIF`と`pymatgen.Structure`の両方で構造データを取得します。
60 - `pymatgen.SpacegroupAnalyzer`を使用して構造を対称化し、空間群情報、格子パラメータ、等価サイトを取得します。
61 - 対称化された構造情報（格子、空間群、原子サイト）を基に、新しいCIFファイルを生成して保存します。
62 - ログファイルへのリダイレクトやエラーハンドリングも行います。
63
64 :returns: None
65 """
66 logfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-out.txt"])
67 print(f"Open logfile [{logfile}]")
68 app.redirect(targets = ["stdout", logfile], mode = 'w')
69
70 convCIFfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-symmetrized.cif"])
71
72 print("")
73 print(f"input : {infile}")
74 print(f"log file: {logfile}")
75 print(f"output symmetrized CIF file: {convCIFfile}")
76
77#=====================================
78# Ready by tkProg
79#=====================================
80 print("")
81 print(f"Read [{infile}] for tkCIF")
82 cif = tkCIF()
83 cifdata = cif.ReadCIF(infile, find_valid_structure = True)
84 if cifdata is None:
85 app.terminate(f"Error: Can not read [{infile}]", pause = True)
86 exit()
87
88 cifdata.Print()
89 cry_source = cifdata.GetCrystal()
90# cry_source.print_inf()
91
92#=====================================
93# Ready by pymatgen
94#=====================================
95 print("")
96 print(f"Read [{infile}] for pymatgen")
97 structure = Structure.from_file(infile)
98 print("Structure:", structure)
99
100 print("")
101 print(f"Symmetrized:")
102 analyzer = SpacegroupAnalyzer(structure)
103 symmetrized_structure = analyzer.get_symmetrized_structure()
104 symmetrized_structure_dict = symmetrized_structure.as_dict()
105
106 spg_name, spg_num = symmetrized_structure.get_space_group_info()
107 symmetry_ops = analyzer.get_symmetry_operations()
108 symmetry_matrix = [op.as_dict()['matrix'] for op in symmetry_ops]
109 nsym = len(symmetry_matrix)
110 print("Space group: ", spg_name, spg_num)
111 print(f"Number of symmetry options: {nsym}")
112# print("Symmetry operatoins: ", symmetry_matrices)
113
114 print("")
115 print("Lattice: ")
116 lattice_inf = symmetrized_structure.lattice
117 a, b, c, alpha, beta, gamma = lattice_inf.parameters
118 volume = lattice_inf.volume
119 aij = lattice_inf.matrix
120 print(" Lattice parameters:", a, b, c, alpha, beta, gamma)
121 print(" Matrix: ", aij)
122 print(" Volume: ", volume)
123
124 print("")
125 print("Equivalent sites:")
126 eq_sites = symmetrized_structure.equivalent_sites
127 for sites in eq_sites:
128 composition = sites[0].species
129 for atom_name in composition.keys():
130 print(" ", atom_name, sites[0].frac_coords, " occ=", composition[atom_name])
131# print("dir(sites[0])=", dir(sites[0]))
132# for site in sites:
133# print(" ", site)
134
135 print("")
136 print("All sites:")
137 sites_inf = symmetrized_structure.sites
138 for site in sites_inf:
139 composition = site.species
140 for atom_name in composition.keys():
141 print(" ", atom_name, site.frac_coords, " occ=", composition[atom_name])
142
143#=====================================
144# Save
145#=====================================
146 print("")
147 print(f"Save the symmetrized structure to [{convCIFfile}]")
148# This does not save the symmetrized structure. So implement with tkCrytal
149# symmetrized_structure.to(filename = convCIFfile, fmt = "cif")
150 out = tkFile(convCIFfile, 'w')
151 if out.fp is None:
152 app.terminate(f"Error: Can not write to [{convCIFfile}]", pause = True)
153
154 out.write(f"data_{cry_source.path}\n")
155 out.write(f"_cell_length_a {a}\n")
156 out.write(f"_cell_length_b {b}\n")
157 out.write(f"_cell_length_c {c}\n")
158 out.write(f"_cell_angle_alpha {alpha}\n")
159 out.write(f"_cell_angle_beta {beta}\n")
160 out.write(f"_cell_angle_gamma {gamma}\n")
161 out.write(f"_cell_volume {volume}\n")
162 out.write(f"\n")
163 out.write(f"_symmetry_space_group_name_H-M '{spg_name}'\n")
164 out.write(f"_symmetry_Int_Tables_number {spg_num}\n")
165 out.write(f"_chemical_formula_structural {cry_source.ChemicalFormula}\n")
166 out.write(f"_chemical_formula_sum {cry_source.ChemicalFormula}\n")
167 out.write(f"_cell_formula_units_Z {cry_source.ChemicalFormulaUnit}\n")
168
169 def to_str(v):
170 """
171 浮動小数点数をCIF形式の分数の文字列に変換する。
172
173 詳細説明:
174 - 特定の浮動小数点値（0, ±1, ±0.5, ±0.25, ±1/3, ±2/3, ±1/6, ±5/6など）を、それぞれの分数表現（例: '1/2', '-1/3'）に変換します。
175 - 正の整数または単純な分数ではない場合は、入力値を文字列として返します。
176 - 0.0の場合はNoneを返します。
177
178 :param v: float
179 変換する浮動小数点数。
180 :returns: Optional[str]
181 変換された文字列、または変換できない場合は元の値、または0.0の場合はNone。
182 """
183 eps = 1.0e-6
184 if v == 0.0:
185 return None
186 if v == 1.0:
187 return '+'
188 if v == -1.0:
189 return '-'
190 if abs(v - 0.5) < eps:
191 return '1/2'
192 if abs(v + 0.5) < eps:
193 return '-1/2'
194 if abs(v - 0.25) < eps:
195 return '1/4'
196 if abs(v + 0.25) < eps:
197 return '-1/4'
198 if abs(v - 1.0/3.0) < eps:
199 return '1/3'
200 if abs(v + 1.0/3.0) < eps:
201 return '-1/3'
202 if abs(v - 2.0/3.0) < eps:
203 return '2/3'
204 if abs(v + 2.0/3.0) < eps:
205 return '-2/3'
206 if abs(v - 1.0/6.0) < eps:
207 return '1/6'
208 if abs(v + 1.0/6.0) < eps:
209 return '-1/6'
210 if abs(v - 5.0/6.0) < eps:
211 return '5/6'
212 if abs(v + 5.0/6.0) < eps:
213 return '-5/6'
214 if v > 0.0:
215 return re.sub(r'^\+', '', v)
216
217 return v
218
219 def vector_to_str(v):
220 """
221 3DベクトルをCIF形式の対称操作文字列（例: 'x', 'x+y', '-x+z'）に変換する。
222
223 詳細説明:
224 - 入力された3要素ベクトル（[x, y, z]）の各要素を`to_str`関数で変換し、'x', 'y', 'z'を連結して対称操作の文字列を生成します。
225 - 例えば、[1, 0, 0] は 'x' に、[0, 1, 0] は 'y' に、[0, 0, -1] は '-z' に、[1, 1, 0] は 'x+y' になります。
226 - 先頭の'+'は除去されます。
227
228 :param v: List[float]
229 3要素の浮動小数点数ベクトル。
230 :returns: str
231 対称操作を表す文字列。
232 """
233 sx = to_str(v[0])
234 sy = to_str(v[1])
235 sz = to_str(v[2])
236 s = ''
237 if sx:
238 s += sx + 'x'
239 if sy:
240 s += sy + 'y'
241 if sz:
242 s += sz + 'z'
243
244 return re.sub(r'^\+', '', s)
245
246 def matrix_to_str(m):
247 """
248 3x3行列をCIF形式の対称操作文字列のタプル（'x', 'y', 'z'成分）に変換する。
249
250 詳細説明:
251 - 入力された3x3行列の各行（ベクトル）を`vector_to_str`関数に渡し、それぞれの変換結果をタプルとして返します。
252 - これは、CIFの `_symmetry_equiv_pos_as_xyz` エントリの 'x, y, z' 部分を生成するために使用されます。
253
254 :param m: List[List[float]]
255 3x3の浮動小数点数行列。
256 :returns: Tuple[str, str, str]
257 'x', 'y', 'z'成分を表す文字列のタプル。
258 """
259 x = vector_to_str(m[0])
260 y = vector_to_str(m[1])
261 z = vector_to_str(m[2])
262
263 return x, y, z
264
265 out.write(f"\n")
266 out.write(f"loop_\n")
267 out.write(f" _symmetry_equiv_pos_site_id\n")
268 out.write(f" _symmetry_equiv_pos_as_xyz\n")
269 for i in range(nsym):
270 x, y, z = matrix_to_str(symmetry_matrix[i])
271 out.write(f" {i+1:2} '{x}, {y}, {z}'\n")
272
273 out.write(f"\n")
274 out.write(f"loop_\n")
275 out.write(f" _atom_site_type_symbol\n")
276# out.write(f" _atom_site_label\n")
277 out.write(f" _atom_site_fract_x\n")
278 out.write(f" _atom_site_fract_y\n")
279 out.write(f" _atom_site_fract_z\n")
280 out.write(f" _atom_site_occupancy\n")
281 for sites in eq_sites:
282 composition = sites[0].species
283 for atom_name in composition.keys():
284 pos = sites[0].frac_coords
285 occ = composition[atom_name]
286
287 out.write(f" {atom_name}")
288 out.write(f" {pos[0]:12.8f} {pos[1]:12.8f} {pos[2]:12.8f}")
289 out.write(f" {occ}\n")
290
291 out.close()
292
293
294 app.terminate(pause = True)
295
296
297if __name__ == "__main__":
298 main()