"""
tklib連携構造解析プログラム
概要:
このスクリプトは、tklibライブラリ（tkCIF, tkCrystal, tkAtomType）とPymatgenライブラリを連携させ、
CIFファイルの読み込み、構造情報の表示、密度計算、結合価数和（BVS）計算、
動径分布関数（RDF）/累積配位数（RCN）計算、Ewald法によるマデルングポテンシャル計算、
構造因子計算、Pymatgenによる対称化構造のCIF出力など、一連の構造解析機能を提供します。
詳細説明:
コマンドライン引数としてCIFファイルパスを受け取り、指定されたCIFファイルに対して様々な解析機能を実施します。
各機能ブロックは独立しており、構造解析のステップバイステップの実行例として設計されています。
tklibとPymatgenの機能比較や連携のデモンストレーションとしても利用できます。
関連リンク:
:doc:`template_tklib_usage` (もし関連する外部ドキュメントがあればここに記述)
"""
import argparse
import sys
import numpy as np
from pprint import pprint
from typing import Tuple, Optional, List, Dict, Any
import re
import math
import os
from numpy import sin, cos, tan, pi, exp, log, log10, sqrt
from scipy.special import erfc
from itertools import product
# --- tklibのインポート ---
# IMPORTANT: tklibのディレクトリがsys.pathに含まれている必要があります
# 例: sys.path.append("d:/git/tkProg/tklib/python")
try:
from tklib.tkfile import tkFile
from tklib.tkcrystal.tkcif import tkCIF
from tklib.tkcrystal.tkcrystal import tkCrystal
from tklib.tkcrystal.tkatomtype import tkAtomType
from tklib.tkcrystal.tkcif2pymatgen import tkcrystal_to_pmg_structure
except ImportError:
print("エラー: tklibライブラリが見つかりません。", file=sys.stderr)
print("sys.pathにtklibのパスが正しく追加されているか確認してください。", file=sys.stderr)
sys.exit(1)
# --- pymatgenのインポート ---
try:
from pymatgen.core import Structure
from pymatgen.core.lattice import Lattice
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
except ImportError:
print("エラー: pymatgenライブラリが見つかりません。", file=sys.stderr)
sys.exit(1)
# --- グローバル定数 (Ewald法 & 構造因子用) ---
# Ewald法に必要な物理定数
e_charge = 1.60218e-19 # C
e0 = 8.854418782e-12 # F/m
# Keはクーロンポテンシャルの係数 e / (4 * pi * e0)
Ke = e_charge / (4.0 * pi * e0)
# X-ray波長 (nm)
wavelength = {
'MoKa1': 0.070926 # nm
, 'MoKa2': 0.071354
, 'MoKa' : 0.071069
, 'CuKa1': 0.15405
, 'CuKa2': 0.15443
, 'CuKa' : 0.15418
, 'CuKb1': 0.13810
, 'CuKb2': 0.13922
, 'CuKb' : 0.13847
}
# --- ヘルパー関数 (共通) ---
[ドキュメント]
def print_matrix(m):
"""
3x3行列を見やすいフォーマットで出力します。
Parameters:
:param m: np.ndarray: 出力する3x3のNumPy行列。
Returns:
:returns: None
"""
print(f"| {m[0][0]:8.4f} {m[0][1]:8.4f} {m[0][2]:8.4f} |")
print(f"| {m[1][0]:8.4f} {m[1][1]:8.4f} {m[1][2]:8.4f} |")
print(f"| {m[2][0]:8.4f} {m[2][1]:8.4f} {m[2][2]:8.4f} |")
# --- BVS計算ヘルパー関数 (機能 5用) ---
[ドキュメント]
def read_bv_params(path: str) -> Optional[Dict[str, Any]]:
"""
指定されたパスからBVSパラメータファイルを読み込み、辞書として返します。
詳細説明:
ファイルが存在しない場合はエラーメッセージを表示しNoneを返します。
ファイル内の '_valence_param_details' セクション以降の行をパースし、
原子ペア (例: 'Atom1:Atom2') をキーとして 'r0', 'B', 'ref' のパラメータを辞書に格納します。
行が '#' で始まる場合はセクションの終わりとみなします。
Parameters:
:param path: str: BVSパラメータファイルへのパス。
Returns:
:returns: Optional[Dict[str, Any]]: BVSパラメータの辞書、またはファイルが見つからない場合にNone。
"""
bv_params = {}
if not os.path.exists(path):
print(f"エラー: BVSパラメータファイルが見つかりません: [{path}]", file=sys.stderr)
return None
print(f"BVSパラメータファイルを読み込み中: [{os.path.basename(path)}]")
with open(path, 'r', encoding='utf-8') as f:
lines = f.readlines()
reading = False
for line in lines:
line = line.strip()
if not line:
continue
if line == '_valence_param_details':
reading = True
continue
if reading:
if line.startswith('#'):
break
parts = line.split()
if len(parts) < 7:
continue
atom1, z1_str, atom2, z2_str, r0_str, b_str, ref = parts[:7]
try:
r0 = float(r0_str); B = float(b_str)
except ValueError:
continue
key = f"{atom1}:{atom2}"
if key not in bv_params:
bv_params[key] = {
'r0': r0, 'B': B, 'ref': ref
}
return bv_params
# --- Pymatgen構造に基づくヘルパー関数 (機能 4用) ---
[ドキュメント]
def site_effective_occupancy_and_mass_amu(site) -> Tuple[float, float]:
"""
Pymatgen Siteの有効占有数と総質量を計算します。
詳細説明:
サイト内の各原子種の占有数と原子質量を合算し、サイト全体の有効占有数と総質量（amu単位）を計算します。
Parameters:
:param site: pymatgen.core.Site: 計算対象のPymatgen Siteオブジェクト。
Returns:
:returns: Tuple[float, float]: 有効占有数と総質量 (amu単位) のタプル。
"""
occ = 0.0
mass = 0.0
for sp, amt in site.species.items():
occ += float(amt)
amu = getattr(sp, "atomic_mass", 0.0)
mass += float(amt) * float(amu if amu is not None else 0.0)
return occ, mass
[ドキュメント]
def structure_effective_counts(struct: Structure) -> Tuple[float, float, int]:
"""
Pymatgen Structureの有効原子数、総質量、サイト数を計算します。
詳細説明:
構造内の全てのサイトについて、`site_effective_occupancy_and_mass_amu` を用いて
有効占有数と総質量を合算し、さらにサイトの総数も返します。
Parameters:
:param struct: Structure: 計算対象のPymatgen Structureオブジェクト。
Returns:
:returns: Tuple[float, float, int]: 有効原子数、総質量 (amu単位)、サイト数のタプル。
"""
eff_atoms = 0.0
total_mass = 0.0
for site in struct.sites:
occ, mass = site_effective_occupancy_and_mass_amu(site)
eff_atoms += occ
total_mass += mass
return eff_atoms, total_mass, len(struct.sites)
[ドキュメント]
def report_pmg_structure_density(struct: Structure, label: str) -> Tuple[float, float]:
"""
Pymatgen Structureの密度を計算し、その結果をレポートとして表示します。
詳細説明:
構造の体積、有効原子数、総質量（amuからgに変換）を計算し、
原子密度（atoms/ų）と質量密度（g/cm³）を算出して出力します。
Parameters:
:param struct: Structure: 計算対象のPymatgen Structureオブジェクト。
:param label: str: 出力メッセージのラベルとして使用される文字列。
Returns:
:returns: Tuple[float, float]: 原子密度 (atoms/ų) と質量密度 (g/cm³) のタプル。
"""
a, b, c = struct.lattice.abc
volume = struct.volume
eff_atoms, total_mass_amu, n_sites = structure_effective_counts(struct)
total_mass_g = total_mass_amu * 1.66053906660e-24
volume_cm3 = volume * 1.0e-24
rho_atom = eff_atoms / volume if volume > 0 else float("nan") # atoms/Å^3
rho_mass = total_mass_g / volume_cm3 if volume_cm3 > 0 else float("nan") # g/cm^3
print(f"\n[{label} Structure (Pymatgen)]")
print(f" 体積 (Volume) : {volume:.6f} ų")
print(f" 有効原子数 (Sigma occ) : {eff_atoms:.6f}")
print(f" 総質量 (Mass) : {total_mass_amu:.6f} amu")
print(f" 原子密度 (Atomic rho): {rho_atom:.6f} atoms/ų")
print(f" 質量密度 (Mass rho): {rho_mass:.6f} g/cm³")
return rho_atom, rho_mass
[ドキュメント]
def change_basis_preserving_geometry(structure: Structure,
T: np.ndarray,
t: Optional[np.ndarray] = None,
xyz_tol: float = 1e-6) -> Structure:
"""
Pymatgen Structureオブジェクトの基底を変換し、重複サイトをマージします。
詳細説明:
指定された変換行列 `T` と並進ベクトル `t` を用いて、格子ベクトルと原子の分率座標を変換します。
変換後の座標は周期境界条件を考慮して `[0, 1)` の範囲に正規化されます。
その後、`xyz_tol` の許容誤差内で同一と見なされるサイトはマージされ、新しいStructureオブジェクトが生成されます。
Parameters:
:param structure: Structure: 基底変換を行う元のPymatgen Structureオブジェクト。
:param T: np.ndarray: 3x3の変換行列。新しい基底ベクトルを定義します。
:param t: Optional[np.ndarray]: 変換後の座標に適用する並進ベクトル (1x3)。デフォルトはNone (並進なし)。
:param xyz_tol: float: サイトの重複を判定するための座標許容誤差。
Returns:
:returns: Structure: 変換され、重複がマージされた新しいPymatgen Structureオブジェクト。
Raises:
:raises ValueError: 変換行列 `T` が特異である (行列式がゼロに近い) 場合。
"""
print(" => Transformation matrix [T]:")
print_matrix(T)
detT = float(np.linalg.det(T))
if abs(detT) < 1e-12:
raise ValueError("Transformation matrix is singular (det ~ 0).")
V = np.array(structure.lattice.matrix, dtype=float)
Vp = T @ V
Tinv = np.linalg.inv(T)
frac = np.array([site.frac_coords for site in structure.sites], dtype=float)
frac_p = frac @ Tinv
if t is not None:
frac_p = frac_p + t.reshape(1,3)
frac_p = np.mod(frac_p, 1.0)
# 同一位置の原子のうち１つだけ残してkept_species/coordsに残す
round_ndp = max(0, int(-np.log10(max(xyz_tol, 1e-12)))) + 1
kept_coords: List[np.ndarray] = []
kept_species = []
seen = set()
for fc, site in zip(frac_p, structure.sites):
key = (round(fc[0]+xyz_tol, round_ndp), round(fc[1]+xyz_tol, round_ndp), round(fc[2]+xyz_tol, round_ndp),
site.species)
if key in seen: continue
seen.add(key)
kept_species.append(site.species)
kept_coords.append(fc)
new_lattice = Lattice(Vp)
new_structure = Structure(new_lattice, kept_species, kept_coords, coords_are_cartesian=False)
new_structure.sort()
removed = len(structure.sites) - len(new_structure.sites)
if removed > 0:
print(f" [Info] Merged {removed} duplicate sites due to det(T)={detT:.6f}.")
return new_structure
# --- 解析機能ブロック 1: tkCIFによる読み込みと基本情報 ---
[ドキュメント]
def analyze_tkcif_basic(filepath: str) -> Optional[tkCrystal]:
"""
tkCIFを使用してCIFファイルを読み込み、基本情報を表示し、tkCrystalオブジェクトを返します。
詳細説明:
指定されたパスのCIFファイルを `tkCIF.ReadCIF` メソッドで読み込みます。
読み込みに成功した場合、`cifdata.Print()` でCIFデータの内容をコンソールに出力し、
そのデータから `tkCrystal` オブジェクトを取得して返します。
ファイル読み込みや `tkCrystal` オブジェクト取得に失敗した場合はエラーメッセージを表示しNoneを返します。
Parameters:
:param filepath: str: 読み込むCIFファイルへのパス。
Returns:
:returns: Optional[tkCrystal]: 読み込まれた構造を表すtkCrystalオブジェクト、
またはエラーが発生した場合にNone。
"""
print("\n" + "="*50)
print("--- 1. tkCIFによるファイル読み込みとCIFデータ ---")
print("="*50)
cif = tkCIF()
cif.debug = 0
try:
cifdata = cif.ReadCIF(filepath, find_valid_structure=1)
except Exception as e:
print(f"エラー: tkCIF.ReadCIFの実行中にエラーが発生しました: {e}", file=sys.stderr)
return None
print(f"ファイル: [{filepath}]")
print("\n[tkCIFData.Print()出力]")
cifdata.Print()
cry = cifdata.GetCrystal()
if cry is None:
print("\nエラー: tkCIFDataからtkCrystalオブジェクトを取得できませんでした。", file=sys.stderr)
return None
return cry
# --- 解析機能ブロック 2: tkCrystalによる構造情報と密度計算 ---
# --- LaTeX修正: 単位から \text{} を削除 ---
[ドキュメント]
def analyze_tkcrystal_structure(cry: tkCrystal):
"""
tkCrystalオブジェクトを用いて構造情報（格子、サイト）と密度を計算し表示します。
詳細説明:
`tkCrystal.PrintInf()` を呼び出して構造の概要情報を出力します。
その後、`cry.Density()` と `cry.AtomDensity()` を用いて質量密度と原子密度を計算し、
その結果を表示するとともに、一般的な密度の範囲内にあるかどうかの簡単な検証を行います。
Parameters:
:param cry: tkCrystal: 解析対象のtkCrystalオブジェクト。
Returns:
:returns: None
"""
print("\n" + "="*50)
print("--- 2. tkCrystalによる構造情報と密度 ---")
print("="*50)
# 1. 構造情報の表示 (PrintInf())
print("\n[tkCrystal.PrintInf()出力]")
print("==============================================")
cry.PrintInf() # 構造の概要（格子、対称操作、サイトなど）を出力
print("==============================================")
# 2. 密度の計算と検証
d = cry.Density() # g/cm3
ad = cry.AtomDensity() # /cm3
print(f"\n[密度計算結果]")
print(f" 質量密度 (Density): {d:.6f} g/cm³") # LaTeX修正
print(f" 原子密度 (AtomDensity): {ad:.2e} atoms/cm³") # LaTeX修正
# 3. 密度の検証 (元のプログラムのロジックを再現)
if not (0.8 < d < 20.0):
print(f"\n❌ 検証失敗: 質量密度が異常な範囲 ({d:.6f} g/cm³) です。")
if not (1.0e22 < ad < 10.0e23):
print(f"\n❌ 検証失敗: 原子密度が異常な範囲 ({ad:.2e} /cm³) です。")
# --- 解析機能ブロック 3: tkAtomTypeの利用（テスト） ---
[ドキュメント]
def test_tkatomtype():
"""
tkAtomTypeを使用して原子の情報を取得するテストです。
詳細説明:
`tkAtomType` クラスのインスタンスを作成し、特定の元素（ここでは'Au'）の原子情報を
`GetAtomInformation` メソッドで取得し、その質量や半径を表示します。
エラーが発生した場合はメッセージを出力します。
Parameters:
:param : なし
Returns:
:returns: None
"""
print("\n" + "="*50)
print("--- 3. tkAtomTypeによる原子情報テスト ---")
print("="*50)
atom = tkAtomType()
test_element = 'Au'
try:
inf = atom.GetAtomInformation(test_element)
print(f"[{test_element}] の原子情報:")
print(f" Atomic Mass: {inf.get('AtomicMass', 'N/A')}")
print(f" Covalent Radius: {inf.get('CovRadius', 'N/A')}")
print(f" Ion Radius (VI-Coord): {inf.get('IonRadiusVI', 'N/A')}")
except Exception as e:
print(f"エラー: tkAtomTypeの実行中にエラーが発生しました: {e}", file=sys.stderr)
# --- 解析機能ブロック 4: Pymatgen連携による単位格子変換と密度検証 ---
# --- LaTeX修正: 単位から \text{} を削除 ---
# --- 解析機能ブロック 5: BVS（結合価数和）計算 ---
# --- LaTeX修正: 単位から \text{} を削除 ---
[ドキュメント]
def calculate_bond_valence_sums(cry: tkCrystal, max_r: float = 2.7):
"""
tkCrystalのサイト情報と自作のBVSパラメータを用いて結合価数和を計算します。
詳細説明:
環境変数 'tkProg_Root' からBVSパラメータファイルのパスを構築し、ファイルを読み込みます。
読み込まれたパラメータと結晶構造情報（非対称サイト、展開サイト）を使用して、
各非対称サイトからの最短原子間距離と結合価数（Bond Valence）を計算し、
その合計である結合価数和（Bond Valence Sum, BVS）を出力します。
結合距離が `max_r` を超える結合は無視されます。
Parameters:
:param cry: tkCrystal: BVS計算の対象となるtkCrystalオブジェクト。
:param max_r: float: 結合と見なす最大原子間距離（Å単位）。
Returns:
:returns: None
"""
print("\n" + "="*50)
print("--- 5. BVS（結合価数和）計算 ---")
print("="*50)
# --- 1. 環境変数とデータベースパスのチェック ---
tkProg_Root = os.environ.get('tkProg_Root')
if not tkProg_Root:
print("❌ エラー: 環境変数 'tkProg_Root' が設定されていません。", file=sys.stderr)
print(" BVS計算をスキップします。", file=sys.stderr)
return
db_path = os.path.join(tkProg_Root, 'tkdb', 'Databases', 'bvparm2020.cif')
# --- 2. BVSパラメータの読み込み ---
bv_params = read_bv_params(db_path)
if bv_params is None:
return
# --- 3. BVS計算の実行 ---
# 非対称単位の代表サイトを取得
asym_sites = cry.AtomSiteList()
expanded_sites = cry.ExpandedAtomSiteList()
print(f"\nBVS計算を開始します (カットオフ距離: {max_r:.2f} A)")
for i, site1 in enumerate(asym_sites):
name1 = site1.AtomNameOnly()
charge1 = site1.Charge()
occ1 = site1.Occupancy()
pos1 = site1.Position(1) # [0, 1) 範囲の分率座標
bvs = 0.0
# サイト名と期待される電荷の表示
charge_str = f"({charge1:.0f}+)" if charge1 > 0 else f"({charge1:.0f}-)" if charge1 < 0 else ""
print(f"\n{i:2d}: {name1}{charge_str} (Occ: {occ1:.2f})")
for site2 in expanded_sites:
name2 = site2.AtomNameOnly()
occ2 = site2.Occupancy()
pos2 = site2.Position(1)
# 2サイト間の最短原子間距離を計算 (AllowZero=Falseで同一サイトを無視)
d = cry.GetNearestInterAtomicDistance(pos1, pos2, AllowZero=False, irange=[1, 1, 1])
if d is None or d > max_r:
continue
# BVSパラメータの取得 (順方向または逆方向)
fwd = f"{name1}:{name2}"
bwd = f"{name2}:{name1}"
params = bv_params.get(fwd) or bv_params.get(bwd)
if not params:
print(f" ⚠️ Warning: BV params for {fwd} or {bwd} not found, skipping bond at {d:.3f} A.")
continue
r0 = params['r0']; B = params['B']
if B == 0: continue
# 結合価数 (Bond Valence) の計算
bv = np.exp((r0 - d) / B)
# サイト2の占有率 $occ_2$ を乗算して加算
bvs += bv * occ2
# 詳細出力
print(f" → {name2:<2} r0={r0:.3f}, B={B:.2f}, d={d:.3f}, bv={bv:.3f} ({params['ref']})")
# 結合価数和の最終出力
print(f"\n ◆ Bond Valence Sum (BVS): {bvs:.3f}\n")
# --- 解析機能ブロック 6: RDF/RCN 計算とデータ出力 ---
[ドキュメント]
def make_matrix_of_size(dims: List[int], dtype=float, defval=0.0):
"""
指定された次元の行列（リストのリスト）を作成するヘルパー関数。
numpy.zerosの代替として、ネストされたPythonリストで初期化された行列を生成します。
Parameters:
:param dims: List[int]: 行列の次元を指定する整数のリスト (例: [dim1, dim2, ...])。
:param dtype: type: 要素のデータ型（この実装では直接は使用されませんが、Docstringには含めます）。
:param defval: Any: 各要素の初期値。
Returns:
:returns: Any: 指定された次元と初期値で作成された行列（リストのリスト）。
"""
if not dims:
return defval
if len(dims) == 1:
return [defval] * dims[0]
return [make_matrix_of_size(dims[1:], dtype, defval) for _ in range(dims[0])]
[ドキュメント]
def CalRDFi(cry: tkCrystal, Rmin: float, Rmax: float, nR: int, atom0, iSite: int, normalize: bool = True):
"""
単一の非対称サイト atom0 から見た全原子に対するタイプ別RDF（動径分布関数）を計算します。
詳細説明:
指定された中心原子 `atom0` から、結晶内の全ての展開された原子サイトまでの距離を計算します。
周期境界条件を考慮し、距離範囲 `Rmin` から `Rmax` を `nR` 個のビンに分割し、
各ビンの原子数（占有率を考慮）を原子タイプ別に集計します。
`normalize` がTrueの場合、結果は距離ビン幅で規格化され、RDF(r) = dN/dr となります。
Parameters:
:param cry: tkCrystal: 計算対象のtkCrystalオブジェクト。
:param Rmin: float: 距離計算の最小値（Å単位）。
:param Rmax: float: 距離計算の最大値（Å単位）。
:param nR: int: 距離ビンの数。
:param atom0: tkAtomSite: 中心となる非対称サイトオブジェクト。
:param iSite: int: `atom0` の非対称サイトリストにおけるインデックス。
:param normalize: bool: 結果を距離ビン幅で規格化するかどうか。Trueの場合、RDF(r)として機能します。
Returns:
:returns: Tuple[List[float], List[List[float]]]:
`xR` (距離配列、各ビンの中心距離のリスト) と
`RDF` (タイプ別のRDFヒストグラム、`RDF[iAtomType1][iR]` の形式)。
"""
print(f"Calculating RDF for site #{iSite}: ", end = '')
Rstep = Rmax / (nR - 1)
xR = [Rstep * k for k in range(nR)]
# PBC探索範囲（Rmaxに基づいて決定）
nLatticeX, nLatticeY, nLatticeZ = cry.GetLatticeRange(Rmax)
AtomTypes = cry.AtomTypeList()
nTypes = len(AtomTypes)
AllSites = cry.ExpandedAtomSiteList()
nSites = len(AllSites)
name0 = atom0.AtomNameOnly()
iAtomType0 = atom0.iAtomType()
x0, y0, z0 = atom0.Position(1)
occ = atom0.Occupancy()
print(f"{name0:<4}[{iAtomType0}] ({x0:8.4f}, {y0:8.4f}, {z0:8.4f}) occ={occ:8.4f}")
# RDF[iAtomType1][iR]
RDF = make_matrix_of_size([nTypes, nR], defval=0.0)
# 全ての展開されたサイトに対して距離を計算
for atom1 in AllSites:
iAtomType1 = atom1.iAtomType()
x1, y1, z1 = atom1.Position(1)
occ1 = atom1.Occupancy()
# 周期境界条件 (PBC) を考慮した探索
for iz in range(-nLatticeZ, nLatticeZ+1):
for iy in range(-nLatticeY, nLatticeY+1):
for ix in range(-nLatticeX, nLatticeX+1):
# サイト0 (中心) と、サイト1の周期像 [x1+ix, y1+iy, z1+iz] の距離
dis = cry.GetInterAtomicDistance([x0, y0, z0], [x1+ix, y1+iy, z1+iz])
# Rmin (自明な距離) を除く & Rmax 以下
if dis < Rmin or Rmax < dis:
continue
idx = int(dis / Rstep)
# ビンに占有率を加算 (Rstepで割らない)
if idx < nR:
RDF[iAtomType1][idx] += occ1
if normalize:
# 距離ビン幅で規格化して RDF(r) = dN/dr にする
for itype in range(nTypes):
for iR in range(nR):
RDF[itype][iR] /= Rstep
return xR, RDF # xR (距離配列) と RDF (タイプ別ヒストグラム)
[ドキュメント]
def CalculateRDFs(cry: tkCrystal, Rmin: float, Rmax: float, nR: int):
"""
全非対称サイトに対するRDFを計算し、それを原子種ベースで平均化してRDF/RCNを求めます。
詳細説明:
まず、`CalRDFi` を使用して、結晶内の全ての非対称サイトをそれぞれ中心としたタイプ別RDFを個別に計算します。
次に、同じ原子種タイプに属する非対称サイトのRDFを平均化し、原子種別平均RDF (`RDFs`) を算出します。
最後に、`RDFs` を累積積分することで原子種別累積配位数 (`RCNs`) を計算します。
Parameters:
:param cry: tkCrystal: 計算対象のtkCrystalオブジェクト。
:param Rmin: float: 距離計算の最小値（Å単位）。
:param Rmax: float: 距離計算の最大値（Å単位）。
:param nR: int: 距離ビンの数。
Returns:
:returns: Tuple[List[float], np.ndarray, np.ndarray]:
`xR` (距離配列)、
`RDFs` (原子種別平均RDF、`RDFs[iType][iR]` の形式)、
`RCNs` (原子種別RCN、`RCNs[iType][iR]` の形式)。
"""
print("\n[ステップ 1] 全非対称サイトに対する個別RDFを計算")
AtomTypes = cry.AtomTypeList()
nTypes = len(AtomTypes)
AsymSites = cry.AtomSiteList()
nAsymSites = len(AsymSites)
# RDF_site[iSite][iAtomType1][iR]
RDF_site = make_matrix_of_size([nAsymSites])
xR = []
for isite in range(nAsymSites):
atom0 = AsymSites[isite]
xR, RDF0 = CalRDFi(cry, Rmin, Rmax, nR, atom0, isite, normalize=True)
RDF_site[isite] = RDF0
print(f"\n[ステップ 2] 原子種別平均RDF (RDFs) の計算")
# RDFs[iType][iR]
RDFs = np.zeros([nTypes, nR], dtype=float)
mult = np.zeros(nTypes, dtype=int)
# 非対称サイトごとに、属する原子種タイプに加算
for isite in range(nAsymSites):
atom1 = AsymSites[isite]
itype1 = atom1.iAtomType()
mult[itype1] += 1 # 各原子種タイプの非対称サイト数をカウント
for itype in range(nTypes):
# サイト isite の中心から見た、タイプ itype の原子のRDFを積算
RDFs[itype1] += np.array(RDF_site[isite][itype])
# 平均化
print(" 原子種別マルチプリシティ (非対称サイト数):")
for itype in range(nTypes):
type_name = AtomTypes[itype].AtomType()
print(f" {itype}: {type_name}: {mult[itype]}")
if mult[itype] > 0:
RDFs[itype] /= mult[itype]
# RCNs[iType][iR] の計算 (累積積分)
RCNs = np.zeros([nTypes, nR], dtype=float)
Rstep = xR[1] - xR[0]
for itype in range(nTypes):
# 累積積分 (Rstepを掛けて $\int \frac{dN}{dR} dR$ を計算)
RCNs[itype] = np.cumsum(RDFs[itype]) * Rstep
return xR, RDFs, RCNs
[ドキュメント]
def calculate_cn_rdfs(cry: tkCrystal, Rmin: float, Rmax: int, nR: int, nCN: int):
"""
CN（配位数）に基づいた累積配位数分布（RCN）とRDFを計算します。
詳細説明:
まず、各非対称サイトからの総原子数（raw count）のRDFとRCNを計算します。
次に、各サイトのRCNに基づき、特定の配位数 `nCN` 以上となる距離範囲で、
原子タイプ別に累積されたRDF (`CNRDFs`) とRCN (`CNRCNs`) を計算します。
`CNRDFs` はCNが `icn` 以上となる距離ビンにおける原子数、`CNRCNs` はその累積積分です。
Parameters:
:param cry: tkCrystal: 計算対象のtkCrystalオブジェクト。
:param Rmin: float: 距離計算の最小値（Å単位）。
:param Rmax: float: 距離計算の最大値（Å単位）。
:param nR: int: 距離ビンの数。
:param nCN: int: 計算する最大配位数。`iCN` は 1 から `nCN` まで使用されます。
Returns:
:returns: Tuple[List[float], List[List[List[float]]], List[List[List[float]]]]:
`xR_cn` (距離配列)、
`CNRDFs` (CNベースのRDF、`CNRDFs[iType][iCN][iR]` の形式)、
`CNRCNs` (CNベースのRCN、`CNRCNs[iType][iCN][iR]` の形式)。
"""
print(f"\n[ステップ 3] CN (配位数) ベースのRCN/RDF（CN {nCN} or greater相当）を計算")
AtomTypes = cry.AtomTypeList()
nTypes = len(AtomTypes)
AsymSites = cry.AtomSiteList()
nAsymSites = len(AsymSites)
# RDFs_total[iSite][iR]: サイト iSite から見た全原子のRDF (Rstepで割ってない生カウント)
RDFs_raw_total = np.zeros([nAsymSites, nR], dtype=float)
for isite in range(nAsymSites):
atom = AsymSites[isite]
# normalize=False で生カウントを取得
xR, RDF0_types = CalRDFi(cry, Rmin, Rmax, nR, atom, isite, normalize=False)
for itype in range(nTypes):
RDFs_raw_total[isite] += np.array(RDF0_types[itype])
# RCNs_total[iSite][iR]: 累積配位数 (生カウントを累積)
RCNs_total = np.zeros([nAsymSites, nR], dtype=float)
for isite in range(nAsymSites):
RCNs_total[isite] = np.cumsum(RDFs_raw_total[isite])
# CNRDFs[iType][iCN][iR] & CNRCNs[iType][iCN][iR]
# iCN は 1から nCN まで使用するため、サイズは nCN+1
CNRDFs = make_matrix_of_size([nTypes, nCN + 1, nR], defval=0.0)
CNRCNs = make_matrix_of_size([nTypes, nCN + 1, nR], defval=0.0)
Rstep = xR[1] - xR[0]
for isite in range(nAsymSites):
atom = AsymSites[isite]
itype = atom.iAtomType()
mult = atom.Multiplicity() * atom.Occupancy() # 原子種の重み
icn = 1
for iR in range(nR):
# iR のビンまでの累積配位数が icn 以上であれば、そのビンは CN>=icn に寄与する
if RCNs_total[isite][iR] >= icn:
# CNRDFs は、CNが icn 以上の寄与を積算（原子種タイプ iType の重み mult を乗算）
# RDFs_raw_total[isite][iR] は、このRでのサイト isite からの原子数 (dN)
CNRDFs[itype][icn][iR] += RDFs_raw_total[isite][iR] * mult
# 次のCNをチェック
icn += 1
if icn > nCN:
break
# CNRCNs の計算 (CNRDFs の累積積分)
for itype in range(nTypes):
for icn in range(1, nCN + 1):
# CNRDFs は既に生カウントのまま累積されているため、Rstepで割ってRDF(r)に変換し、積分する
CNRDF_normalized = np.array(CNRDFs[itype][icn]) / Rstep
CNRCNs[itype][icn] = np.cumsum(CNRDF_normalized) * Rstep
return xR, CNRDFs, CNRCNs
# --- 解析機能ブロック 7: Ewald法によるマデルングポテンシャル計算 ---
[ドキュメント]
def calculate_madelung_potential(cry: tkCrystal, ew_alpha: float = 0.3, prec: float = 1.0e-5):
"""
Ewald法を用いて、単位格子内の最初の非対称サイトにおけるマデルングポテンシャルを計算します。
詳細説明:
Ewald法は、クーロン相互作用の実空間項、逆空間項、自己項の3つの部分に分けて計算することで、
無限に広がる結晶の静電ポテンシャルを収束させます。
この関数では、ターゲットとなる最初の非対称サイトに対してこれらの項を計算し、合計してマデルングポテンシャル
（Joule単位およびeV単位）とマデルング定数を出力します。
実空間および逆空間のカットオフは、指定された精度 `prec` に基づいて決定されます。
Parameters:
:param cry: tkCrystal: 計算対象のtkCrystalオブジェクト。
:param ew_alpha: float: Ewaldパラメータ $\alpha$。実空間と逆空間の寄与のバランスを制御します。
:param prec: float: 計算の目標精度。結果の相対誤差に影響します。
Returns:
:returns: None
"""
print("\n" + "="*50)
print("--- 7. Ewald法によるマデルングポテンシャル計算 ---")
print("="*50)
sites = cry.AtomSiteList()
if not sites:
print("❌ エラー: 構造内に非対称サイトがありません。", file=sys.stderr)
return
# ターゲットサイト (最初の非対称サイト)
target_site = sites[0]
pos_i = np.array(target_site.Position(1)) # 分率座標
qi = target_site.Charge() * target_site.Occupancy() # 有効電荷
# 格子情報
# 修正: gij, Rgij = cry.Metrics() を使用
gij, Rgij = cry.Metrics()
volume = cry.Volume() * 1.0e-30 # Å^3 -> m^3
Rvolume = cry.ReciprocalVolume() # Å^-3
print(f"ターゲットサイト: {target_site.AtomName()} (電荷: {qi:.4f})")
print(f"Ewaldパラメータ: alpha={ew_alpha}, precision={prec}")
print(f"単位セル体積: {cry.Volume():.4g} ų")
if abs(qi) < 1e-6:
print("警告: ターゲットサイトの有効電荷がゼロに近いため、ポテンシャルはゼロに近くなります。", file=sys.stderr)
return
# --- 1. 実空間項 (UC1) の計算 ---
# リアル空間のカットオフ距離 RDmax の決定
norder = -log10(prec)
rdmax = (2.26 + 0.26 * norder) / ew_alpha
print(f" Real space cut-off: RDmax = {rdmax:.4f} A")
# 探索範囲 $\mathbf{n}=(n_x, n_y, n_z)$ の最大値決定
# $\mathbf{n}$ の最大値を決定するため、格子ベクトル長と $R_{Dmax}$ を使用
nrmax = [int(rdmax / cry.LatticeParameters()[i]) + 2 for i in range(3)]
UC1 = 0.0
# 実空間和 (UC1) の計算
for pos_j_frac in cry.AtomSiteList():
qj = pos_j_frac.Charge() * pos_j_frac.Occupancy()
pos_j = np.array(pos_j_frac.Position(1)) # 分率座標
# 格子点の探索
for nz in range(-nrmax[2], nrmax[2]+1):
for ny in range(-nrmax[1], nrmax[1]+1):
for nx in range(-nrmax[0], nrmax[0]+1):
n_vec = np.array([nx, ny, nz])
# サイト間距離の計算
rij = cry.GetInterAtomicDistance(pos_i, pos_j + n_vec, AllowZero=True, irange=[0, 0, 0])
# サイト $i$ 自身（$r_{ii}=0$）をスキップ (自己項はUC3で扱う)
is_self = (abs(rij) < 1.0e-12) and (qj == qi) and np.allclose(pos_i, pos_j) and np.allclose(n_vec, [0, 0, 0])
if is_self:
continue
if rij > rdmax:
continue
# 実空間項: $V_1 = \sum_{j,\mathbf{n} \ne (i,\mathbf{0})} q_j \frac{\mathrm{erfc}(\alpha r_{ij})}{r_{ij}}$
erfcar = erfc(ew_alpha * rij)
UC1 += qj * erfcar / (rij * 1.0e-10)
# --- 2. 逆空間項 (UC2) の計算 ---
G2max = ew_alpha**2 / pi**2 * (-log(prec))
print(f" Reciprocal space cut-off: G2max = {G2max:.4f} A⁻²")
# 逆空間の探索範囲 $\mathbf{h}=(h, k, l)$ の最大値決定
# Rgijは対角要素ではないため、正確には対角要素のみで近似するのは不正確だが、元のコードの意図を尊重
hmax = [int(sqrt(G2max / Rgij[i][i])) + 2 for i in range(3)]
UC2 = 0.0
Krec = 1.0 / (pi * volume)
Kexp = pi * pi / ew_alpha / ew_alpha
# 逆格子点 $\mathbf{h}$ の探索
for l in range(-hmax[2], hmax[2]+1):
for k in range(-hmax[1], hmax[1]+1):
for h in range(-hmax[0], hmax[0]+1):
h_vec = np.array([h, k, l])
# $G^2 = |\mathbf{G}|^2 = \mathbf{h} \cdot G \cdot \mathbf{h}^T$
G2 = np.dot(h_vec, np.dot(Rgij, h_vec))
if G2 < 1.0e-12 or G2 > G2max: # $\mathbf{h}=(0,0,0)$ とカットオフを超えるものをスキップ
continue
# 構造因子 $S(\mathbf{G})$ の実数部 $\text{Re}[ S(\mathbf{G}) e^{-i \mathbf{G} \cdot \mathbf{r}_i} ]$
cossum_j = 0.0
sinsum_j = 0.0
for site_j in cry.AtomSiteList(): # 全てのサイト $j$
qj = site_j.Charge() * site_j.Occupancy()
pos_j = np.array(site_j.Position(1))
# $\phi_j = 2\pi \mathbf{h} \cdot \mathbf{r}_j$
phi_j = 2.0 * pi * np.dot(h_vec, pos_j)
cossum_j += qj * cos(phi_j)
sinsum_j += qj * sin(phi_j)
phi_i = 2.0 * pi * np.dot(h_vec, pos_i)
cosphi_i = cos(phi_i)
sinphi_i = sin(phi_i)
# $\sum_j q_j \cos(\mathbf{G} \cdot (\mathbf{r}_i - \mathbf{r}_j))$
fcal = cosphi_i * cossum_j + sinphi_i * sinsum_j
# $e^{-\pi^2 G^2 / \alpha^2}$ / $G^2$
expg = exp(-Kexp * G2) / (G2 * 1.0e+20)
# $\mathbf{h}$ と $-\mathbf{h}$ の項を一度に計算するため $2.0$ 倍
UC2 += Krec * expg * fcal * 2.0
# --- 3. 自己項 (UC3) の計算 ---
# 自己項: $V_3 = -q_i \frac{2\alpha}{\sqrt{\pi}}$
UC3 = qi * 2.0 * (ew_alpha * 1.0e10) / sqrt(pi)
# --- 4. ポテンシャルの合計 ---
MP = UC1 + UC2 - UC3
# --- 5. 結果の表示 ---
print("\n[Ewald法 結果]")
print(f" 実空間和 (UC1) : {UC1:.6g} V")
print(f" 逆空間和 (UC2) : {UC2:.6g} V")
print(f" 自己項 (UC3) : {UC3:.6g} V")
# Madelung Potential (J)
mad_pot_J = Ke * MP
print(f" マデルングポテンシャル: {mad_pot_J:.6g} J")
# Madelung Potential (eV)
mad_pot_eV = Ke / e_charge * MP
print(f" マデルングポテンシャル: {mad_pot_eV:.6g} eV")
# Madelung Constant (NaCl $\alpha \approx 1.7475$)
a_ang = cry.LatticeParameters()[0] # 最初の格子定数 a (Å)
mad_const = 0.5 * MP / abs(qi) * (a_ang * 1.0e-10) * e0 * (4 * pi)
print(f" マデルング定数: {mad_const:.8g} (格子定数 a で規格化)")
# --- 解析機能ブロック 8: 構造因子（Fhkl）と原子散乱因子の計算 ---
[ドキュメント]
def calculate_structure_factor(cry: tkCrystal, xray_source_key: str = 'CuKa1', hkl_test: List[int] = [1, 0, 0]):
"""
指定されたX線波長と$hkl$での構造因子 $|F_{hkl}|$ と原子散乱因子を計算・表示します。
詳細説明:
まず、X線源キーから波長を取得し、指定されたミラー指数 (hkl) と波長に基づいて
面間隔 $d_{hkl}$、散乱ベクトル $s = \sin\theta/\lambda$、および回折角 $2\theta$ を計算します。
次に、結晶内の各原子タイプについて、散乱ベクトル $s$ に応じた原子散乱因子 $f(s)$ を `tkAtomType.asf()` で求め、
最後に `tkCrystal.Fhkl()` を用いて結晶構造因子 $F_{hkl}$ の複素数とその絶対値を計算し、出力します。
Parameters:
:param cry: tkCrystal: 構造因子計算の対象となるtkCrystalオブジェクト。
:param xray_source_key: str: 使用するX線源のキー（例: 'CuKa1'）。`wavelength` 辞書に定義されている必要があります。
:param hkl_test: List[int]: 計算するミラー指数 [h, k, l] のリスト。
Returns:
:returns: None
"""
print("\n" + "="*50)
print("--- 8. 構造因子（Fhkl）と原子散乱因子の計算 ---")
print("="*50)
wl_nm = wavelength.get(xray_source_key)
if wl_nm is None:
print(f"❌ エラー: X線源 '{xray_source_key}' の波長データが見つかりません。", file=sys.stderr)
return
# tkCrystalは距離を Å (オングストローム) で扱うため、wl_nmを wl_A (Å) に変換
wl_A = wl_nm * 10.0
print(f"X線源: {xray_source_key} (波長: {wl_nm:.5f} nm = {wl_A:.5f} A)")
print(f"ターゲットミラー指数 (hkl): ({hkl_test[0]}, {hkl_test[1]}, {hkl_test[2]})")
# $d_{hkl}$, $s = \sin\theta/\lambda$, $2\theta$ の計算
try:
# cry.DiffractionAngleは $\text{wl}$ を Å、 $s$ を $\text{A}^{-1}$、 $Q2$ を $\text{rad}$ で返す前提
dhkl_A, s_A_inv, two_theta_rad = cry.DiffractionAngle(wl_A, *hkl_test)
except Exception as e:
print(f"❌ エラー: 回折角の計算に失敗しました: {e}", file=sys.stderr)
return
if dhkl_A is None:
print("警告: 該当するhklでの回折は理論上起こりません (dhklが計算不能)。")
return
s_nm_inv = s_A_inv / 10.0 # $s$ を $\text{nm}^{-1}$ に変換 (tkAtomType.asf の入力単位に合わせる)
two_theta_deg = two_theta_rad * 180.0 / pi
print("\n[回折情報]")
print(f" 面間隔 $d_{{hkl}}$: {dhkl_A:.6g} A")
print(f" 散乱ベクトル $s$: {s_nm_inv:.6g} nm⁻¹")
print(f" 回折角 $2\\theta$: {two_theta_deg:.6g} deg")
# 1. 原子散乱因子 $f(s)$ の計算
print("\n[原子散乱因子 $f(s)$]")
AtomTypes = cry.AtomTypeList()
for at in AtomTypes:
atom_name = at.AtomTypeOnly()
# $s$ は $\text{nm}^{-1}$ 単位で渡す
f_s = at.asf(s_nm_inv) # f(s)
# f_prime, f_double_prime = at.AnomalousScatteringFactor(wl_nm) # $f'$, $f''$
f_prime, f_double_prime = 0.0, 0.0
# 複合散乱因子 $f = f(s) + f' + i f''$
f_complex = f_s + f_prime + 1j * f_double_prime
print(f" {atom_name:<4} -> $f(s)={f_s:.4f}$, $f'={f_prime:.4f}$, $f''={f_double_prime:.4f}$")
print(f" 複合 $f$: {f_complex.real:.6g} + j{f_complex.imag:.6g}")
# 2. 構造因子 $F_{hkl}$ の計算
# cry.Fhkl は $s$ を $\text{A}^{-1}$ 単位で渡す前提の可能性が高いため、 $s\_A\_inv$ を使用
try:
Fhkl_complex = cry.Fhkl(s_A_inv, *hkl_test)
# Fhklの絶対値と位相
Fhkl_abs = np.abs(Fhkl_complex)
print(f"\n[結晶構造因子 $F_{{hkl}}({hkl_test[0]}, {hkl_test[1]}, {hkl_test[2]})]$")
print(f" $F_{{hkl}}$ (複素数): {Fhkl_complex.real:.6g} + j{Fhkl_complex.imag:.6g}")
print(f" $|F_{{hkl}}|$ (絶対値): {Fhkl_abs:.6g}")
except Exception as e:
print(f"❌ エラー: 構造因子の計算に失敗しました: {e}", file=sys.stderr)
# --- 解析機能ブロック 9: Pymatgen対称化後の結果をtkFileでCIF出力 ---
[ドキュメント]
def output_symmetrized_cif(cif_filepath: str, cry_source: tkCrystal, sym_tol: float = 1.0e-3):
"""
tkCrystalをPymatgenで対称化し、その結果をtkFileを使ってCIF形式で保存します。
詳細説明:
`tkCrystal` オブジェクトをPymatgenの `Structure` に変換し、`SpacegroupAnalyzer` を用いて
対称化された構造 (`symmetrized_structure`) を取得します。
対称化された構造から空間群情報、格子パラメータ、対称操作、非等価サイト情報を抽出し、
新しいCIFファイルに整形して書き出します。対称操作のXYZ表現の生成にはヘルパー関数が使用されます。
Parameters:
:param cif_filepath: str: 元のCIFファイルのパス。出力ファイル名の生成に使用されます。
:param cry_source: tkCrystal: 対称化の対象となるtkCrystalオブジェクト。
:param sym_tol: float: PymatgenのSpacegroupAnalyzerが対称性を検出する際の許容誤差（Å単位）。
Returns:
:returns: Optional[str]: 出力されたCIFファイルのパス、またはエラーが発生した場合にNone。
"""
print("\n" + "="*50)
print("--- 9. Pymatgen対称化後の結果をtkFileでCIF出力 ---")
print("="*50)
# 1. Pymatgen Structureへの変換と対称化
try:
structure = tkcrystal_to_pmg_structure(cry_source)
analyzer = SpacegroupAnalyzer(structure, symprec=sym_tol)
symmetrized_structure = analyzer.get_symmetrized_structure()
spg_name, spg_num = symmetrized_structure.get_space_group_info()
lattice = symmetrized_structure.lattice
eq_sites = symmetrized_structure.equivalent_sites
# SymmOpの行列を取得（CIF出力用）
symmetry_ops = analyzer.get_symmetry_operations()
symmetry_matrix = [op.rotation_matrix for op in symmetry_ops]
symmetry_translation = [op.translation_vector for op in symmetry_ops]
nsym = len(symmetry_ops)
print(f"✅ 対称化に成功。空間群: {spg_name} #{spg_num}")
print(f" 対称操作数: {nsym}")
except Exception as e:
print(f"❌ エラー: 対称化またはPymatgen処理に失敗しました: {e}", file=sys.stderr)
return
# 出力ファイル名
header, ext = os.path.splitext(cif_filepath)
convCIFfile = f"{header}_symmetrized.cif"
print(f"\n[出力ファイル]: {convCIFfile}")
# 2. CIF整形ヘルパー関数 (元のプログラムから抽出)
eps = 1.0e-6
def to_str(v):
"""
数値を分数または特殊な文字列表現に変換します。
Parameters:
:param v: float: 変換対象の数値。
Returns:
:returns: Optional[str]: 変換された文字列 ('+', '-', '1/2' など)、または数値文字列。
数値がeps以下の場合None。
"""
if abs(v) < eps: return None
if abs(v - 1.0) < eps: return '+'
if abs(v + 1.0) < eps: return '-'
if abs(v - 0.5) < eps: return '1/2'
if abs(v + 0.5) < eps: return '-1/2'
if abs(v - 1.0/3.0) < eps: return '1/3'
if abs(v + 1.0/3.0) < eps: return '-1/3'
if abs(v - 2.0/3.0) < eps: return '2/3'
if abs(v + 2.0/3.0) < eps: return '-2/3'
# その他の場合は数値として返す
return f"{v:0.8f}"
def matrix_to_xyz_str(matrix_row, translation_component):
"""
回転行列の行と並進成分からCIF形式の 'x, y, z' 形式の文字列を生成します。
Parameters:
:param matrix_row: np.ndarray: 回転行列の1行 (例: [1, 0, 0])。
:param translation_component: float: 対応する並進成分。
Returns:
:returns: str: CIFの `_symmetry_equiv_pos_as_xyz` に記述される形式の文字列。
"""
labels = ['x', 'y', 'z']
terms = []
for i, coeff in enumerate(matrix_row):
if abs(coeff) > eps:
coeff_str = to_str(coeff)
if coeff_str is None: continue # 0.0はスキップ
if coeff_str in ('+', '-'):
term = f"{coeff_str}{labels[i]}"
elif coeff_str.startswith('-'):
term = f"{coeff_str[1:]}{labels[i]}"
else:
term = f"{coeff_str}{labels[i]}"
# 最初でなければ '+' をつける（負符号は既に付いている）
if terms and not term.startswith('-'):
terms.append(f"+{term}")
else:
terms.append(term)
# 並進成分の追加
if abs(translation_component) > eps:
t_str = to_str(translation_component)
if t_str is None: t_str = ''
if terms and t_str not in ('+', '-') and not t_str.startswith('-'):
terms.append(f"+{t_str}")
else:
terms.append(t_str)
return "".join(terms).replace('++', '+').replace('+-', '-').lstrip('+').lstrip('-') or '0'
# 3. tkFileでのCIF書き出し
try:
out = tkFile(convCIFfile, 'w')
except Exception as e:
print(f"❌ エラー: ファイル [{convCIFfile}] の書き込みに失敗しました: {e}", file=sys.stderr)
return
# ヘッダー情報の書き出し
a, b, c, alpha, beta, gamma = lattice.parameters
out.write(f"data_{os.path.splitext(os.path.basename(cif_filepath))[0]}_symmetrized\n")
out.write(f"_cell_length_a {a:.8f}\n")
out.write(f"_cell_length_b {b:.8f}\n")
out.write(f"_cell_length_c {c:.8f}\n")
out.write(f"_cell_angle_alpha {alpha:.4f}\n")
out.write(f"_cell_angle_beta {beta:.4f}\n")
out.write(f"_cell_angle_gamma {gamma:.4f}\n")
out.write(f"_cell_volume {lattice.volume:.8f}\n")
out.write(f"\n")
out.write(f"_symmetry_space_group_name_H-M '{spg_name}'\n")
out.write(f"_symmetry_Int_Tables_number {spg_num}\n")
out.write(f"_chemical_formula_structural {cry_source.ChemicalFormula}\n")
out.write(f"_chemical_formula_sum {cry_source.ChemicalFormula}\n")
out.write(f"_cell_formula_units_Z {cry_source.ChemicalFormulaUnit}\n")
# 対称操作の書き出し
out.write(f"\n")
out.write(f"loop_\n")
out.write(f" _symmetry_equiv_pos_site_id\n")
out.write(f" _symmetry_equiv_pos_as_xyz\n")
for i in range(nsym):
rot_row = symmetry_matrix[i]
trans_comp = symmetry_translation[i]
x_eq = matrix_to_xyz_str(rot_row[0], trans_comp[0])
y_eq = matrix_to_xyz_str(rot_row[1], trans_comp[1])
z_eq = matrix_to_xyz_str(rot_row[2], trans_comp[2])
out.write(f" {i+1:2} '{x_eq}, {y_eq}, {z_eq}'\n")
# サイト情報の書き出し (非等価サイト)
out.write(f"\n")
out.write(f"loop_\n")
out.write(f" _atom_site_type_symbol\n")
out.write(f" _atom_site_fract_x\n")
out.write(f" _atom_site_fract_y\n")
out.write(f" _atom_site_fract_z\n")
out.write(f" _atom_site_occupancy\n")
for sites in eq_sites:
# 等価サイト群の最初の要素（代表サイト）を使用
rep_site = sites[0]
pos = rep_site.frac_coords
# 複合原子種の処理
for specie, occ in rep_site.species.items():
atom_name = specie.symbol
out.write(f" {atom_name:<5}")
out.write(f" {pos[0]:12.8f} {pos[1]:12.8f} {pos[2]:12.8f}")
out.write(f" {occ:.8f}\n")
out.close()
print("✅ 対称化されたCIFを正常に保存しました。")
return convCIFfile
# --- メイン関数 ---
[ドキュメント]
def main():
"""
CLI引数からCIFファイルを受け取り、tkcif, tkcrystalの機能を順次実行します。
詳細説明:
コマンドライン引数としてCIFファイルのパスを受け取ります。
引数が指定されない場合は、デフォルトの 'ZnO.cif' を使用します。
その後、以下の構造解析機能ブロックを順番に実行します。
1. tkCIFによるCIF読み込みと基本情報表示
2. tkCrystalによる構造情報と密度計算
3. tkAtomTypeによる原子情報テスト
4. Pymatgen連携による単位格子変換と密度検証
5. 結合価数和（BVS）計算
6. 動径分布関数（RDF）/累積配位数（RCN）計算
7. Ewald法によるマデルングポテンシャル計算
8. 構造因子（Fhkl）と原子散乱因子の計算
9. Pymatgen対称化後の結果をtkFileでCIF出力
Parameters:
:param : なし
Returns:
:returns: None
"""
default_cif = "ZnO.cif"
parser = argparse.ArgumentParser(
description="tklib (tkCIF/tkCrystal)の機能を実行する構造解析プログラム。",
formatter_class=argparse.RawTextHelpFormatter
)
parser.add_argument(
"cif_file",
type=str,
nargs='?',
default=default_cif,
help=f"解析するCIFファイルのパス。\n(デフォルト: {default_cif} を想定)"
)
args = parser.parse_args()
print(f"tklib参照プログラムを開始します。対象ファイル: {args.cif_file}\n")
# --- 機能ブロック 1: CIF読み込みとCrystalオブジェクト取得 ---
cry_object = analyze_tkcif_basic(args.cif_file)
if cry_object:
# --- 機能ブロック 2: Crystalオブジェクトの解析と密度 ---
analyze_tkcrystal_structure(cry_object)
# --- 機能ブロック 4: Pymatgen連携による単位格子変換 ---
perform_pmg_conversion_on_tkcrystal(cry_object, sym_tol=0.01, density_eps=1.0e-4)
# --- 機能ブロック 5: BVS計算 ---
calculate_bond_valence_sums(cry_object, max_r=2.7)
# --- 機能ブロック 6: RDF/RCN 計算とデータ出力 ---
print("\n" + "="*50)
print("--- 6. RDF/RCN 計算とデータ出力 ---")
print("="*50)
Rmin, Rmax, nR = 0.1, 5.0, 1001
nCN = 6
# 6-1. 原子種ベースのRDF/RCN計算
print("\n[テスト 1] 原子種ベースの平均RDF/RCNの計算 (Atom mode相当)")
try:
xR, RDFs, RCNs = CalculateRDFs(cry_object, Rmin, Rmax, nR)
print(f"✅ RDFs/RCNsの計算に成功しました。Rmax={Rmax:.1f} A, nR={nR}")
# 結果の概要表示（ZnOのZnサイトの最初のピークをチェック）
zn_type_index = cry_object.AtomTypeList()[0].iAtomType() # 最初のタイプをZnと仮定
# RDFsが空でないことを確認
if RDFs.size > 0 and np.any(RDFs[zn_type_index] > 0):
# RCNsがリストのリストであることを確認
if isinstance(RCNs, np.ndarray) and RCNs.ndim == 2:
first_rcn_peak = RCNs[zn_type_index][np.argmax(RDFs[zn_type_index] > 0)]
print(f" ZnタイプのRDFの最初の非ゼロRCN: {first_rcn_peak:.4f}")
else:
print(" RCNsの形状が予期されたものではありませんでした。")
else:
print(" RDFデータが空またはゼロです。")
except Exception as e:
print(f"❌ RDFs/RCNsの計算中にエラーが発生しました: {e}", file=sys.stderr)
# 6-2. CNベースのRDF/RCN計算
print(f"\n[テスト 2] CN (配位数) ベースのCNRDF/CNRCNの計算 (CN {nCN} or greater相当)")
try:
xR_cn, CNRDFs, CNRCNs = calculate_cn_rdfs(cry_object, Rmin, Rmax, nR, nCN)
print(f"✅ CNRDFs/CNRCNsの計算に成功しました。Max CN={nCN}")
# 結果の概要表示（ZnサイトのCN>=4の累積数をチェック）
zn_type_index = cry_object.AtomTypeList()[0].iAtomType()
# CNRCNsが3次元配列であることを確認
if len(CNRCNs) > zn_type_index and len(CNRCNs[zn_type_index]) > 4:
# CNRCNs[itype][icn][iR]
cn4_final = CNRCNs[zn_type_index][4][-1] # Type 0, CN=4, 最後のR
print(f" ZnタイプのCN 4 or greater の最大累積数: {cn4_final:.4f}")
except Exception as e:
print(f"❌ CNRDFs/CNRCNsの計算中にエラーが発生しました: {e}", file=sys.stderr)
# --- 機能ブロック 7: Ewald法によるマデルングポテンシャル計算 ---
# ZnS型構造（ZnO）では電荷が設定されている必要がある
calculate_madelung_potential(cry_object)
# --- 機能ブロック 8: 構造因子計算 ---
calculate_structure_factor(cry_object, xray_source_key='CuKa1', hkl_test=[1, 0, 0])
# --- 機能ブロック 9: 対称化結果のCIF保存 ---
output_symmetrized_cif(args.cif_file, cry_object)
# --- 機能ブロック 3: tkAtomTypeのテスト ---
test_tkatomtype()
print("\n--- 全ての機能の実行を完了しました。 ---")
if __name__ == "__main__":
main()