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行列を見やすいフォーマットで出力します。"""
    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パラメータファイルを読み込み、辞書として返します。"""
    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の有効占有数と総質量を計算します。"""
    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の有効原子数、総質量、サイト数を計算します。"""
    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の密度をレポートします。"""
    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オブジェクトの基底を変換し、重複サイトをマージします。
    """
    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オブジェクトを返します。
    """
    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オブジェクトを用いて構造情報（格子、サイト）と密度を計算し表示します。
    """
    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を使用して原子の情報を取得するテストです。
    """
    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{} を削除 ---
def perform_pmg_conversion_on_tkcrystal(cry: tkCrystal, sym_tol: float = 0.01, density_eps: float = 1.0e-4):
    """
    tkCrystalをPymatgen Structureに変換し、Primitive Cellへの変換と密度検証を実行します。
    """
    print("\n" + "="*50)
    print("--- 4. Pymatgen連携による単位格子変換と密度検証 ---")
    print("="*50)

    # 1. tkCrystal -> Pymatgen Structure への変換
    try:
        s_orig = tkcrystal_to_pmg_structure(cry)
        print("✅ tkCrystalからPymatgen Structureへの変換に成功しました。")
        print(f"  Pymatgen Structure: {s_orig.formula}, Sites: {len(s_orig)}")
    except Exception as e:
        print(f"❌ エラー: Pymatgen Structureへの変換に失敗しました: {e}", file=sys.stderr)
        return

    # 2. 元のPymatgen Structureのレポート (密度を記録)
    print("\n[ステップ 1] 変換前のPymatgen構造情報と密度計算")
    rho_atom_orig, rho_mass_orig = report_pmg_structure_density(s_orig, "Original")

    # 3. Primitive Cellへの変換
    try:
        print("\n[ステップ 2] Primitive Standard Structureへの変換 (SpacegroupAnalyzerを使用)")
        analyzer = SpacegroupAnalyzer(s_orig, symprec=sym_tol)
        s_conv = analyzer.get_primitive_standard_structure()
    except Exception as e:
        print(f"警告: Standard Structureの取得に失敗。原始格子へのフォールバック: {e}", file=sys.stderr)
        s_conv = s_orig.get_primitive_structure()
        
    # 4. 変換後の構造のレポート
    print("\n[ステップ 3] 変換後のPymatgen構造情報と密度計算")
    rho_atom_new, rho_mass_new = report_pmg_structure_density(s_conv, "Converted (Primitive)")

    # 5. 密度の一致確認
    print(f"\n[ステップ 4] 密度の一貫性チェック (Pymatgen -> Pymatgen)") 

    is_atom_rho_consistent = np.isclose(rho_atom_orig, rho_atom_new, rtol=density_eps)
    if is_atom_rho_consistent:
        print(f"✅ 原子密度は一致しています (許容誤差 eps={density_eps})") 
    else:
        print(f"❌ 警告: 原子密度が一致しませんでした。Original: {rho_atom_orig:.6f}, Converted: {rho_atom_new:.6f}")
        
    is_mass_rho_consistent = np.isclose(rho_mass_orig, rho_mass_new, rtol=density_eps)
    if is_mass_rho_consistent:
        print(f"✅ 質量密度は一致しています (許容誤差 eps={density_eps})") 
    else:
        print(f"❌ 警告: 質量密度が一致しませんでした。Original: {rho_mass_orig:.6f}, Converted: {rho_mass_new:.6f}")

    print(f"\n[Info] Primitive Structureの体積は、元の体積の {s_conv.volume / s_orig.volume:.4f} 倍です。")


# --- 解析機能ブロック 5: BVS（結合価数和）計算 ---
# --- LaTeX修正: 単位から \text{} を削除 ---
def calculate_bond_valence_sums(cry: tkCrystal, max_r: float = 2.7):
    """
    tkCrystalのサイト情報と自作のBVSパラメータを用いて結合価数和を計算します。
    """
    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の代替として使用。
    """
    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（タイプ別）を計算します。
    """
    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を求めます。
    """
    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: float, nR: int, nCN: int):
    """
    CN（配位数）に基づいた累積配位数分布（RCN）とRDFを計算します。
    """
    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法を用いて、単位格子内の最初の非対称サイトにおけるマデルングポテンシャルを計算します。
    """
    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}|$ と原子散乱因子を計算・表示します。
    """
    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形式で保存します。
    """
    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):
        """数値を分数または特殊な文字列表現に変換"""
        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):
        """回転行列の行と並進成分から 'x, y, z' 形式の文字列を生成"""
        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の機能を順次実行します。
    """
    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()