Initial commit: add iodide-PS4 sharing analysis script

2026-06-29 13:14:45 +09:00
commit 79d3f8ff26
8 changed files with 2485231 additions and 0 deletions
@@ -0,0 +1,194 @@
+#!/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}' に保存されました。")