lpsi-iodide-environment-analysis/calc_share.py

#!/usr/bin/env python
# ./calc_share_v2.py -i "*.lammpstrj"
import numpy as np
import argparse
import os
import glob
from collections import Counter
import pandas as pd


class LoadData():
    def __init__(self, trjfile):
        self.trjfile = trjfile
        self.loadtrj()

    def loadtrj(self):
        print(f" -> [処理中] {os.path.basename(self.trjfile)} を読み込み中...")
        with open(self.trjfile) as o:
            d = o.read()
        lines = d.splitlines()    
        self.atoms = int(lines[3])

        x_bnd = [float(x) for x in lines[5].split()]
        y_bnd = [float(x) for x in lines[6].split()]
        z_bnd = [float(x) for x in lines[7].split()]
        self.box_size = np.array(
            [x_bnd[1]-x_bnd[0], y_bnd[1]-y_bnd[0], z_bnd[1]-z_bnd[0]])

        data = np.array(" ".join(lines).split())
        self.data = data.reshape(-1, 28+self.atoms*6)
        self.step = self.data.shape[0]
        data = self.data[:, 28:].reshape(self.step, self.atoms, -1)
        
        self.p_data = data[data[:, :, 2] == "P"].reshape(self.step, -1, 6)
        self.p_xyz = self.p_data[:, :, 3:].astype(float)
        
        self.i_data = data[data[:, :, 2] == "I"].reshape(self.step, -1, 6)
        self.i_xyz = self.i_data[:, :, 3:].astype(float)

        self.s_data = data[data[:, :, 2] == "S"].reshape(self.step, -1, 6)
        self.s_xyz = self.s_data[:, :, 3:].astype(float)
        
        return data

    def get_pbc_distance(self, xyz1, xyz2):
        d = xyz1[:, np.newaxis, :] - xyz2[np.newaxis, :, :]
        d -= self.box_size * np.round(d / self.box_size)
        return np.linalg.norm(d, axis=-1)

    def classify_I(self):
        I_S_CUTOFF = 4.7
        total_counts = Counter()
        total_I_count = 0 
        
        for step in range(self.step):
            p_xyz = self.p_xyz[step]
            s_xyz = self.s_xyz[step]
            i_xyz = self.i_xyz[step]
            
            total_I_count += len(i_xyz)
            
            dist_ps = self.get_pbc_distance(s_xyz, p_xyz)
            s_to_p = np.argmin(dist_ps, axis=1)
            dist_is = self.get_pbc_distance(i_xyz, s_xyz)
            
            for i in range(len(i_xyz)):
                close_s_indices = np.where(dist_is[i] <= I_S_CUTOFF)[0]
                associated_ps4 = s_to_p[close_s_indices]
                ps4_counts = Counter(associated_ps4)
                
                edge = 0
                corner = 0
                for p_idx, s_count in ps4_counts.items():
                    if s_count == 2:
                        edge += 1
                    elif s_count == 1:
                        corner += 1
                
                state = (edge, corner)
                total_counts[state] += 1
                
        avg_counts = {}
        for state, count in total_counts.items():
            column_name = f"edge_{state[0]}_corner_{state[1]}"
            if total_I_count > 0:
                avg_counts[column_name] = count / total_I_count
            else:
                avg_counts[column_name] = 0.0
            
        return avg_counts


if __name__ == "__main__":
    par = argparse.ArgumentParser()
    par.add_argument('-i', '--pattern', default="*.lammpstrj",
                     help='input file pattern (e.g., "*.lammpstrj")')
    par.add_argument('-o', '--output', default="summary_table.tex",
                     help='output tex file name')
    args = par.parse_args()

    filepaths = sorted(glob.glob(args.pattern))
    
    if not filepaths:
        print(f"エラー: パターン '{args.pattern}' に一致するファイルが見つかりませんでした。")
        exit()

    print(f"計 {len(filepaths)} 個のファイルを順番に処理します...")

    all_results = []

    for index, filepath in enumerate(filepaths, start=1):
        filename = os.path.splitext(os.path.basename(filepath))[0]
        print(f"[{index}/{len(filepaths)}] ファイル処理中...")
        
        try:
            calc = LoadData(filepath)
            result = calc.classify_I()
            
            result['Composition'] = filename
            all_results.append(result)
            print(f" -> [完了] {filename} の解析が成功しました。")
            
        except Exception as e:
            print(f" -> [エラー] {filepath} の処理中に問題が発生しました: {e}")

    # --- TeX形式の出力処理 ---
    if all_results:
        print("\nTeXコードの生成中...")
        df = pd.DataFrame(all_results)
        
        # 欠損値を0で埋め、Compositionをインデックスに設定
        df = df.fillna(0)
        df.set_index('Composition', inplace=True)
        
        # 列名（状態）をアルファベット・数字順にソート
        state_cols = sorted(df.columns)
        df = df[state_cols]
        
        # 行と列を入れ替える (Transpose)
        df_t = df.T
        df_t.index.name = 'Label'
        df_t = df_t.reset_index() # Labelをデータ列に戻す
        
        # 列数からフォーマット指定子を生成
        num_cols = len(df_t.columns)
        tabular_format = "wl{2cm}" + "wc{2cm}" * (num_cols - 1)
        
        # TeXテキストの構築
        tex_lines = []
        tex_lines.append("\\begin{table}[H]")
        tex_lines.append(f" \\begin{{tabular}}\n {{{tabular_format}}}")
        tex_lines.append(" \\hline")
        
        # ヘッダー行
        headers = [str(h).replace('_', r'\_') for h in df_t.columns]
        tex_lines.append("  " + " & ".join(headers) + " \\\\")
        tex_lines.append(" \\hline")
        
        # データ行
        for _, row in df_t.iterrows():
            row_vals = []
            for i, val in enumerate(row.values):
                if i == 0:
                    # Label列（状態名）のアンダースコアをエスケープ
                    row_vals.append(str(val).replace('_', r'\_'))
                else:
                    # 数値列の処理（百分率化、四捨五入、1%以下の処理）
                    try:
                        num = float(val)
                        pct = num * 100
                        if pct == 0.0:
                            row_vals.append("0.0\\%")
                        elif pct <= 1.0:
                            # 1%以下の場合は数値ではなく「trace（微量）」という英語にする
                            row_vals.append("trace")
                        else:
                            # 小数点1桁で四捨五入
                            row_vals.append(f"{pct:.1f}\\%")
                    except ValueError:
                        row_vals.append(str(val).replace('_', r'\_'))
            tex_lines.append("  " + " & ".join(row_vals) + " \\\\")
            
        tex_lines.append(" \\end{tabular}")
        tex_lines.append("\\end{table}")
        
        tex_output = "\n".join(tex_lines)
        
        # ファイルへの保存
        with open(args.output, "w", encoding="utf-8") as f:
            f.write(tex_output)
            
        print("\n=== TeX Output ===")
        print(tex_output)
        print(f"\n結果は '{args.output}' に保存されました。")