Initinal commit /lps-lii-ps4-iodide-geometry

.gitignore

@@ -0,0 +1,20 @@
+# Ignore everything by default.
+*
+
+# Keep repository documentation and Git settings.
+!.gitignore
+!README.org
+
+# Keep analysis scripts.
+!dump2cube_edge_corner.py
+!cube2mayavi_edge_corner.py
+
+# Keep the selected example trajectory.
+!050Li3PS4-050LiI_thin100.lammpstrj
+
+# Keep the selected example cube files.
+!050Li3PS4-050LiI_PS4_I_All.cube
+!050Li3PS4-050LiI_PS4_I_zero.cube
+!050Li3PS4-050LiI_PS4_I_corner.cube
+!050Li3PS4-050LiI_PS4_I_edge.cube
+!050Li3PS4-050LiI_PS4_I_three.cube

README.org

@@ -0,0 +1,237 @@
+#+TITLE: PS₄³⁻-Centered I⁻ Geometry Analysis in Li₂S–P₂S₅–LiI Glasses
+#+AUTHOR: Minami Sakuma
+#+OPTIONS: toc:2 num:nil
+
+* Overview
+
+This repository contains Python scripts for analyzing the local
+geometrical relationship between I⁻ ions and PS₄³⁻ tetrahedra in
+Li₂S–P₂S₅–LiI glass trajectories.
+
+For each P atom, nearby I⁻ ions are classified according to the number
+of S atoms from the same PS₄³⁻ tetrahedron located within a specified
+I–S cutoff distance.
+
+The resulting I⁻ spatial-density distributions are exported as Gaussian
+cube files and visualized as three-dimensional isosurfaces using Mayavi.
+
+* Analysis Concept
+
+For a given I⁻ ion and a given PS₄³⁻ tetrahedron, the number of nearby
+S atoms is defined as:
+
+#+begin_example
+s_count = number of S atoms within 4.7 Å from I⁻
+#+end_example
+
+The I⁻ configuration is classified as follows:
+
+| Classification | s_count | Geometrical interpretation |
+|---+---:|---|
+| Zero-type | 0 | I⁻ is not close to any S atom of the reference PS₄³⁻ unit. |
+| Corner-type | 1 | I⁻ is located near one S vertex of the PS₄³⁻ tetrahedron. |
+| Edge-type | 2 | I⁻ is located near two S atoms forming one tetrahedral edge. |
+| Three-type | 3 | I⁻ is located near three S atoms of the same PS₄³⁻ unit. |
+
+The terms corner-type and edge-type describe geometrical proximity only.
+They do not imply bond sharing or atom sharing between I⁻ and PS₄³⁻
+units.
+
+* Workflow
+
+#+begin_example
+LAMMPS trajectory
+      |
+      v
+dump2cube_edge_corner.py
+      |
+      +-- PS₄_I_All.cube
+      +-- PS₄_I_zero.cube
+      +-- PS₄_I_corner.cube
+      +-- PS₄_I_edge.cube
+      +-- PS₄_I_three.cube
+      |
+      v
+cube2mayavi_edge_corner.py
+      |
+      v
+Colored 3D isosurface visualization
+#+end_example
+
+* Files
+
+| File                                 | Description                                                                                          |
+|--------------------------------------+------------------------------------------------------------------------------------------------------|
+| =dump2cube_edge_corner.py=           | Reads a LAMMPS trajectory, aligns PS₄³⁻ units, classifies nearby I⁻ ions, and generates cube files. |
+| =cube2mayavi_edge_corner.py=         | Visualizes multiple classified I⁻ cube files using different colors.                                 |
+| =050Li3PS4-050LiI.lammpstrj=         | Example LAMMPS trajectory for a Li₂S–P₂S₅–LiI glass.                                               |
+| =050Li3PS4-050LiI_PS4_I_All.cube=    | Spatial density of the four nearest I⁻ ions around each P atom.                                      |
+| =050Li3PS4-050LiI_PS4_I_zero.cube=   | Spatial density of zero-type I⁻ configurations.                                                      |
+| =050Li3PS4-050LiI_PS4_I_corner.cube= | Spatial density of corner-type I⁻ configurations.                                                    |
+| =050Li3PS4-050LiI_PS4_I_edge.cube=   | Spatial density of edge-type I⁻ configurations.                                                      |
+| =050Li3PS4-050LiI_PS4_I_three.cube=  | Spatial density of three-type I⁻ configurations.                                                     |
+
+* Requirements
+
+** Python packages
+
+The scripts require Python 3 and the following packages:
+
+- NumPy
+- Mayavi
+- VTK
+- Traits
+- PyQt5 or PySide6
+
+Mayavi is generally easiest to install through conda-forge.
+
+#+begin_src bash
+conda create -n ps4-iodide python=3.10
+conda activate ps4-iodide
+conda install -c conda-forge numpy mayavi pyqt
+#+end_src
+
+Alternatively, NumPy can be installed with pip:
+
+#+begin_src bash
+pip install numpy
+#+end_src
+
+* Input Trajectory Format
+
+=dump2cube_edge_corner.py= expects a LAMMPS trajectory containing the
+following atom columns:
+
+#+begin_example
+ITEM: ATOMS id type element xu yu zu mol
+#+end_example
+
+The trajectory must contain at least the following elements:
+
+- Li
+- P
+- S
+- I
+
+The script assumes that P and S atoms belonging to the same PS₄³⁻
+tetrahedron have the same molecule ID (=mol=).
+
+Periodic boundary conditions are applied when calculating relative
+atomic positions.
+
+* Analysis Procedure
+
+For every trajectory frame, the following procedure is performed.
+
+1. Each P atom is selected as the center of a reference PS₄³⁻ tetrahedron.
+2. The nearest I⁻ ion is used to define the orientation of the PS₄³⁻ unit.
+3. The S atom farthest from the nearest I⁻ ion is aligned with the z axis.
+4. A second S atom is used to fix the rotation around the z axis.
+5. The closest I⁻ ions are rotated into the common PS₄³⁻ reference frame.
+6. I⁻ ions are classified using the number of S atoms within 4.7 Å.
+7. The classified I⁻ coordinates are accumulated over all P atoms and
+   trajectory frames.
+8. Three-dimensional histograms are exported as Gaussian cube files.
+
+* Normalization
+
+The density fields for zero-type, corner-type, edge-type, and three-type
+I⁻ configurations are normalized by the total number of I⁻ coordinates
+used for the overall I⁻ distribution.
+
+Therefore, the relative density of each classified field reflects both:
+
+- the spatial distribution of the configuration, and
+- the relative occurrence of that configuration.
+
+* Usage
+
+** Generate Cube Files
+
+#+begin_src bash
+python dump2cube_edge_corner.py \
+  -i 050Li3PS4-050LiI.lammpstrj \
+  -m 160 160 160 \
+  -cut 8
+#+end_src
+
+Arguments:
+
+| Argument | Description |
+|---+---|
+| =-i=, =--trjfile= | Input LAMMPS trajectory file. |
+| =-m=, =--mesh= | Number of grid points along the x, y, and z directions. |
+| =-cut=, =--cutoff= | Spatial cutoff distance around the reference P atom in Å. |
+
+Expected output files:
+
+#+begin_example
+050Li3PS4-050LiI_PS4_I_All.cube
+050Li3PS4-050LiI_PS4_I_zero.cube
+050Li3PS4-050LiI_PS4_I_corner.cube
+050Li3PS4-050LiI_PS4_I_edge.cube
+050Li3PS4-050LiI_PS4_I_three.cube
+#+end_example
+
+** Visualize Classified I⁻ Density Fields
+
+#+begin_src bash
+python cube2mayavi_edge_corner.py \
+  -i 050Li3PS4-050LiI_PS4_I_zero.cube \
+     050Li3PS4-050LiI_PS4_I_corner.cube \
+     050Li3PS4-050LiI_PS4_I_edge.cube \
+     050Li3PS4-050LiI_PS4_I_three.cube \
+  -iso 1.2e-09
+#+end_src
+
+The script assigns colors based on the file name:
+
+| Classification | Color  |
+|----------------+--------|
+| Zero-type      | Blue   |
+| Corner-type    | Green  |
+| Edge-type      | Red    |
+| Three-type     | Yellow |
+
+The reference PS₄³⁻ tetrahedron is shown with:
+
+| Object    | Color  |
+|-----------+--------|
+| P atom    | Purple |
+| S atoms   | Yellow |
+| P–S bonds | Gray   |
+
+* Output
+
+Each Gaussian cube file contains:
+
+- A reference PS₄³⁻ tetrahedron
+  - P atom at the origin
+  - Four S atoms in the aligned coordinate system
+- A three-dimensional spatial-density field of the selected I⁻ category
+
+The cube files can be visualized using:
+
+- Mayavi
+- VMD
+- ParaView
+- PyMOL
+- Other software supporting Gaussian cube files
+
+* Notes and Limitations
+
+- The I–S cutoff for classification is fixed at 4.7 Å in the script.
+- The cube files represent accumulated spatial-density distributions.
+- The current analysis uses the nearest I⁻ ion to define the orientation
+  of each PS₄³⁻ tetrahedron.
+- The script searches the 15 nearest I⁻ ions around each P atom when
+  classifying zero-, corner-, edge-, and three-type configurations.
+- The output filename is generated from the input trajectory name.
+  The current implementation expects an input filename containing =LiI=.
+- The current trajectory parser is primarily intended for orthorhombic
+  simulation cells.
+- The Mayavi visualization requires a GUI-capable Python environment.
+
+* License
+
+This repository is intended for academic and research use.

cube2mayavi_edge_corner.py

@@ -0,0 +1,191 @@
+#!/usr/bin/env python
+# ./cube2mayavi_edge_corner.py -i 050Li3PS4-050LiI_PS4_I_zero.cube 050Li3PS4-050LiI_PS4_I_corner.cube 050Li3PS4-050LiI_PS4_I_edge.cube 050Li3PS4-050LiI_PS4_I_three.cube -iso 1.2e-09
+import numpy as np
+import os
+from mayavi import mlab
+
+ang_borr = 0.5291772109217
+
+
+def read_cube(filename):
+    """
+    .cubeファイルを読み込み、原子座標と3次元ボクセルデータを返す関数
+    """
+    with open(filename, 'r') as f:
+        # ヘッダーの読み飛ばし（最初の2行はコメント）
+        f.readline()
+        f.readline()
+
+        # 3行目: 原子の数 と 原点座標
+        line = f.readline().split()
+        natoms = abs(int(line[0]))  # 原子の数
+        origin = np.array([float(x) for x in line[1:4]])  # 原点
+
+        # 4-6行目: グリッドの分割数(N)とベクトル(X, Y, Z軸)
+        line = f.readline().split()
+        nx, x_vec = int(line[0]), np.array([float(x) for x in line[1:4]])
+        line = f.readline().split()
+        ny, y_vec = int(line[0]), np.array([float(x) for x in line[1:4]])
+        line = f.readline().split()
+        nz, z_vec = int(line[0]), np.array([float(x) for x in line[1:4]])
+
+        # 原子座標の読み込み
+        atoms = []
+        for _ in range(natoms):
+            line = f.readline().split()
+            atoms.append([float(x) for x in line[2:5]])
+        atoms = np.array(atoms)
+
+        # ボクセルデータの読み込み
+        data = []
+        for line in f:
+            data.extend([float(x) for x in line.split()])
+
+    vol_data = np.array(data).reshape(nx, ny, nz)
+
+    return vol_data, atoms, origin, (x_vec, y_vec, z_vec)
+
+
+def draw_radial_scale(max_dist=10.0, step=1.0, axis='x'):
+    """
+    原点から動径方向に目盛り（定規）を描画する
+    """
+    if axis == 'y':
+        vec = np.array([0, 1, 0])
+        tick_dir = np.array([0.2, 0, 0])
+    elif axis == 'z':
+        vec = np.array([0, 0, 1])
+        tick_dir = np.array([0.2, 0, 0])
+    else:  # default x
+        vec = np.array([1, 0, 0])
+        tick_dir = np.array([0, 0, 0.2])
+
+    # 1. メインの直線を引く
+    end_point = vec * max_dist
+    mlab.plot3d([0, end_point[0]], [0, end_point[1]], [0, end_point[2]],
+                tube_radius=0.02, color=(0, 0, 0))
+
+    # 2. 目盛りと数字を配置
+    for r in np.arange(step, max_dist + step, step):
+        pos = vec * r
+        t_start = pos - tick_dir
+        t_end = pos + tick_dir
+
+        mlab.plot3d([t_start[0], t_end[0]],
+                    [t_start[1], t_end[1]],
+                    [t_start[2], t_end[2]],
+                    tube_radius=0.02, color=(0, 0, 0))
+
+        text_pos = pos + tick_dir * 1.5
+        mlab.text3d(text_pos[0], text_pos[1], text_pos[2],
+                    f"{int(r)}", scale=0.25, color=(0, 0, 0))
+
+
+def change_view(azimuth=30, elevation=60):
+    mlab.view(azimuth=azimuth, elevation=elevation)
+
+
+def get_color_from_filename(filename):
+    """
+    ファイル名から Edge, Corner などの判別を行い、RGBカラーを返す
+    """
+    fname_lower = os.path.basename(filename).lower()
+    if 'zero' in fname_lower:
+        return (0.227, 0.373, 0.804)  # 青色
+    elif 'corner' in fname_lower:
+        return (0.000, 0.733, 0.000)  # 緑色
+    elif 'edge' in fname_lower:
+        return (0.850, 0.058, 0.058)  # 赤色
+    elif 'three' in fname_lower:
+        return (0.943, 0.754, 0.000)  # 黄色
+    else:
+        return (0.5, 0.5, 0.5)        # 該当しない場合はグレー
+
+
+def visualize_cubes(filenames, iso_val=None, cutoff=8):
+    """
+    複数のcubeファイルを読み込んで可視化する
+    """
+    mlab.figure(bgcolor=(1, 1, 1), size=(800, 600))
+    atoms_drawn = False
+
+    for filename in filenames:
+        print(f"Reading {filename}...")
+        vol_data, atoms, origin, vectors = read_cube(filename)
+        atoms = atoms * ang_borr
+        origin_ang = origin * ang_borr
+
+        dx = np.linalg.norm(vectors[0]) * ang_borr
+        dy = np.linalg.norm(vectors[1]) * ang_borr
+        dz = np.linalg.norm(vectors[2]) * ang_borr
+
+        # iso_val が未指定の場合はデータから自動計算
+        current_iso = iso_val
+        if current_iso is None:
+            current_iso = np.std(vol_data) * 2.0
+
+        print(f"  -> Plotting isosurface at value: +/- {current_iso:.4e}")
+
+        # --- カットオフによるデータの球状マスク処理 ---
+        if cutoff is not None:
+            nx, ny, nz = vol_data.shape
+            x = origin_ang[0] + np.arange(nx) * dx
+            y = origin_ang[1] + np.arange(ny) * dy
+            z = origin_ang[2] + np.arange(nz) * dz
+            X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
+            R = np.sqrt(X**2 + Y**2 + Z**2)
+            vol_data[R > cutoff] = 0.0
+
+        print(f"  -> max_value: {np.max(vol_data):.4e}, min_value: {np.min(vol_data):.4e}")
+
+        # --- 確率密度 ---
+        src_in = mlab.pipeline.scalar_field(vol_data)
+        src_in.spacing = [dx, dy, dz]
+        src_in.origin = origin_ang
+
+        # ファイル名から色を取得
+        color = get_color_from_filename(filename)
+
+        # 等値面の描画
+        mlab.pipeline.iso_surface(src_in, contours=[current_iso], opacity=0.3, color=color)
+
+        # --- 原子の描画（最初の1回のみ実行） ---
+        if not atoms_drawn:
+            mlab.points3d(atoms[0, 0], atoms[0, 1], atoms[0, 2],
+                          scale_factor=0.8, color=(0.576, 0.439, 0.8), resolution=20)
+            mlab.points3d(atoms[1:, 0], atoms[1:, 1], atoms[1:, 2],
+                          scale_factor=0.8, color=(1, 1, 0), resolution=20)
+
+            # --- 結合の描画 ---
+            p_coord = atoms[0]
+            s_coords = atoms[1:]
+            for s_coord in s_coords:
+                mlab.plot3d([p_coord[0], s_coord[0]],
+                            [p_coord[1], s_coord[1]],
+                            [p_coord[2], s_coord[2]],
+                            tube_radius=0.1, color=(0.6, 0.6, 0.6), opacity=1.0)
+
+            # 半径の基準となるうすい球の描画
+            mlab.points3d(0, 0, 0, scale_factor=cutoff * 2,
+                          mode='sphere', color=(0, 0, 0), opacity=0.05, resolution=50)
+
+            atoms_drawn = True
+
+
+# --- 実行部分 ---
+if __name__ == "__main__":
+    import argparse
+
+    description = """Visualizes multiple cube files with different colors based on their filenames."""
+    par = argparse.ArgumentParser(description=description)
+
+    # 複数ファイルを受け取れるように nargs="+" を設定（既存のまま）
+    par.add_argument('-i', '--infiles', default=[], required=True, nargs="+",
+                     help='input cube files (e.g., *_I_*.cube)')
+    par.add_argument('-iso', '--iso_val', required=False, type=float,
+                     help='Threshold for isosurface (しきい値)')
+    args = par.parse_args()
+
+    # 複数ファイルをリストとして渡す
+    visualize_cubes(args.infiles, iso_val=args.iso_val)
+    mlab.show()

dump2cube_edge_corner.py

@@ -0,0 +1,333 @@
+#!/usr/bin/env python
+# ./dump2cube_edge_corner.py -i 050Li3PS4-050LiI.lammpstrj -m 160 160 160 -cut 8
+import numpy as np
+import argparse
+import re
+import os
+import time
+
+description = """This is a test program"""
+par = argparse.ArgumentParser(description=description)
+
+par.add_argument('-i', '--trjfile', default="", required=True,
+                 help='input trjfile')
+par.add_argument('-m', '--mesh', default=(30, 30, 30), required=False,
+                 nargs=3, type=int, help='mesh grid')
+par.add_argument('-cut', '--cutoff', default=7, required=False,
+                 type=float, help='mesh grid')
+
+args = par.parse_args()
+ang_borr = 0.5291772109217
+dirname = os.path.dirname(__file__)
+
+
+class LammpsTrj():
+    def __init__(self):
+        """
+        trjファイルを読み込んで、原子数、lattice、ステップ数を格納
+        """
+        with open(args.trjfile) as o:
+            data = o.read().split()
+        self.atoms = int(data[7])  # 原子数
+        data = np.array(data).reshape(-1, self.atoms*7+32)
+        self.lattices = data[:, 17:23].astype(float)  # (step, 6)
+        self.data = data[:, 32:]  # 座標データ
+        self.steps = data.shape[0]
+        self.mesh = args.mesh
+        self.cutoff = args.cutoff
+
+    def setLattice(self, lat):
+        """
+        latticeの形によって操作を分岐
+        self.M, self.M_を作成
+        """
+        if lat.shape[0] == 6:
+            M = np.array([[lat[1]-lat[0], 0, 0],
+                          [0, lat[3]-lat[2], 0],
+                          [0, 0, lat[5]-lat[4]]]).astype(float)
+        if lat.shape[0] == 9:
+            xlo_bound, xhi_bound, xy = lat[0, 0], lat[0, 1], lat[0, 2]
+            ylo_bound, yhi_bound, xz = lat[1, 0], lat[1, 1], lat[1, 2]
+            zlo_bound, zhi_bound, yz = lat[2, 0], lat[2, 1], lat[2, 2]
+            xlo = xlo_bound - np.min([0.0, xy, xz, xy+xz])
+            xhi = xhi_bound - np.max([0.0, xy, xz, xy+xz])
+            ylo = ylo_bound - np.min([0.0, yz])
+            yhi = yhi_bound - np.max([0.0, yz])
+            zlo = zlo_bound
+            zhi = zhi_bound
+            lx = xhi - xlo
+            ly = yhi - ylo
+            lz = zhi - zlo
+            a = lx
+            b = np.sqrt(ly**2 + xy**2)
+            c = np.sqrt(lz**2 + xz**2 + yz**2)
+            alpha = np.arccos((xy*xz + ly*yz)/b/c)
+            beta = np.arccos(xz/c)
+            gamma = np.arccos(xy/b)
+            v1 = [a, 0, 0]
+            v2 = [b*np.cos(gamma), b*np.sin(gamma), 0]
+            v3 = [c*np.cos(beta),
+                  c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/np.sin(gamma),
+                  c*np.sqrt(1+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma)
+                            - np.cos(alpha)**2-np.cos(beta)**2
+                            - np.cos(gamma)**2 / np.sin(gamma))]
+            M = np.array([v1, v2, v3])
+        return M
+
+    def getOnestep(self, step):
+        """
+        引数のstepにおける座標を、原子ごとにself.elems(辞書)に格納する
+        """
+        step_data = self.data[step, :].reshape(self.atoms, -1)  # (atoms, 6)
+        self.M = self.setLattice(self.lattices[step])
+        self.M_ = np.linalg.inv(self.M)
+        self.elems = {}
+        for e in np.unique(step_data[:, 2]):
+            d = step_data[step_data[:, 2] == e, 3:].astype(float)
+            dxyz = d[:, 0:3]
+            dxyz = dxyz @ self.M_
+            dxyz = dxyz - np.floor(dxyz)
+            d[:, 0:3] = dxyz
+            self.elems[e] = d
+        return self.elems
+
+    def makeCube(self, trjfile):
+        """
+        P原子ごとに最も近い4つのI原子を選出し、Edge/Corner分類および回転を行う
+        """
+        self.li_coords_list = []
+        self.i_coords_list = []
+        self.i_zero_coords_list = []
+        self.i_corner_coords_list = []
+        self.i_edge_coords_list = []
+        self.i_three_coords_list = []
+
+        self.li_hist = None
+        self.i_hist = None
+        self.i_hist_zero = None
+        self.i_hist_corner = None
+        self.i_hist_edge = None
+        self.i_hist_three = None
+
+        count_p = 0
+        I_S_CUTOFF = 4.7  # Edge/Corner判定用距離 (Angstrom)
+
+        for step in range(self.steps):
+            print(f"processing {step} step")
+            self.elems = self.getOnestep(step)
+            elems_ = {k: v.copy() for k, v in self.elems.items()}
+
+            for p_data in elems_["P"]:
+
+                # P原子に最も近い4つのI原子を特定
+                i_diffs = elems_["I"][:, 0:3] - p_data[0:3]
+                i_diffs = i_diffs - np.around(i_diffs)
+                i_diffs_abs = i_diffs @ self.M
+                distances = np.linalg.norm(i_diffs_abs, axis=1)
+
+                sorted_indices = np.argsort(distances)
+                closest_Is = elems_["I"][sorted_indices[:4]]
+                closer_Is = elems_["I"][sorted_indices[:15]]
+
+                nearest_I_coords = closest_Is[0]
+
+                # 着目P原子と同じPS4を構成するS原子を取得
+                s_data = elems_["S"][elems_["S"][:, 3] == p_data[3]]
+
+                # --- 距離によるSの配置決定 ---
+                si_diff = s_data[:, 0:3] - nearest_I_coords[0:3]
+                si_diff = si_diff - np.around(si_diff)
+                si_diff = np.linalg.norm(si_diff @ self.M, axis=1)
+
+                n_S_idx = np.argmax(si_diff)
+                s_indices = [0, 1, 2, 3]
+                s_indices.remove(n_S_idx)
+                n2_S_idx = s_indices[0]
+
+                s_xyz = s_data[:, 0:3] - p_data[0:3]
+                s_xyz = s_xyz - np.round(s_xyz)
+                s_xyz = s_xyz @ self.M
+
+                # --- 座標の回転行列を作成 ---
+                theta = np.arctan2(s_xyz[n_S_idx][0], s_xyz[n_S_idx][1])
+                self.Mat_z = np.array([[np.cos(-theta), np.sin(-theta), 0],
+                                       [-np.sin(-theta), np.cos(-theta), 0],
+                                       [0, 0, 1]])
+                s_xyz = (self.Mat_z @ s_xyz.T).T
+
+                theta2 = np.arctan2(s_xyz[n_S_idx][1], s_xyz[n_S_idx][2])
+                self.Mat_x = np.array([[1, 0, 0],
+                                       [0, np.cos(-theta2), np.sin(-theta2)],
+                                       [0, -np.sin(-theta2), np.cos(-theta2)]])
+                s_xyz = (self.Mat_x @ s_xyz.T).T
+
+                theta3 = np.arctan2(s_xyz[n2_S_idx][0], s_xyz[n2_S_idx][1])
+                self.Mat_z2 = np.array([[np.cos(-theta3), np.sin(-theta3), 0],
+                                        [-np.sin(-theta3), np.cos(-theta3), 0],
+                                        [0, 0, 1]])
+                s_xyz = (self.Mat_z2 @ s_xyz.T).T
+                self.s_xyz = s_xyz
+
+                p_xyz = np.array([0, 0, 0])
+                self.ps4_coord = np.vstack((p_xyz, s_xyz)) / ang_borr
+
+                def rotate_coords(coords_fractional):
+                    if len(coords_fractional) == 0:
+                        return np.array([])
+                    rot_xyz = coords_fractional[:, 0:3] - p_data[0:3]
+                    rot_xyz = rot_xyz - np.round(rot_xyz)
+                    rot_xyz = rot_xyz @ self.M
+                    rot_xyz = (self.Mat_z @ rot_xyz.T).T
+                    rot_xyz = (self.Mat_x @ rot_xyz.T).T
+                    rot_xyz = (self.Mat_z2 @ rot_xyz.T).T
+                    return rot_xyz
+
+                # --- I原子のEdge/Corner判定と振り分け ---
+                zero_Is = []
+                edge_Is = []
+                corner_Is = []
+                three_Is = []
+                other_Is = []
+
+                for i_coord in closer_Is:
+                    diff_si = s_data[:, 0:3] - i_coord[0:3]
+                    diff_si = diff_si - np.around(diff_si)
+                    diff_si_abs = diff_si @ self.M
+                    dists = np.linalg.norm(diff_si_abs, axis=1)
+
+                    # 距離がCUTOFF以下のS原子の数をカウント
+                    close_S_count = np.sum(dists <= I_S_CUTOFF)
+
+                    if close_S_count == 0:
+                        zero_Is.append(i_coord)
+                    elif close_S_count == 1:
+                        corner_Is.append(i_coord)
+                    elif close_S_count == 2:
+                        edge_Is.append(i_coord)
+                    elif close_S_count == 3:
+                        three_Is.append(i_coord)
+                    elif close_S_count > 3:
+                        other_Is.append(i_coord)
+
+                # 判定結果に基づいてそれぞれのリストに追加
+                if len(zero_Is) > 0:
+                    self.i_zero_coords_list.extend(
+                        rotate_coords(np.array(zero_Is)).tolist())
+                if len(edge_Is) > 0:
+                    self.i_edge_coords_list.extend(
+                        rotate_coords(np.array(edge_Is)).tolist())
+
+                if len(corner_Is) > 0:
+                    self.i_corner_coords_list.extend(
+                        rotate_coords(np.array(corner_Is)).tolist())
+
+                if len(three_Is) > 0:
+                    self.i_three_coords_list.extend(
+                        rotate_coords(np.array(three_Is)).tolist())
+
+                # 全体のI用リストにも追加
+                self.i_coords_list.extend(rotate_coords(closest_Is).tolist())
+
+                count_p += 1
+
+        print("Calculating Histograms...")
+        volume = (self.cutoff*2)**3
+        bounds = [[-(self.cutoff)/ang_borr, (self.cutoff)/ang_borr]] * 3
+
+        # I (全体のI原子数で割ることで、Edge/Cornerの比率も反映させる)
+        total_I_num = len(self.i_coords_list)
+
+        if total_I_num > 0:
+            # I (All)
+            print("num_I_all: ", total_I_num)
+            i_arr = np.array(self.i_coords_list) / ang_borr
+            i_hist, _ = np.histogramdd(i_arr, bins=self.mesh, range=bounds)
+            self.i_hist = i_hist.ravel() / total_I_num / volume
+
+            # I (zero)
+            if len(self.i_zero_coords_list) > 0:
+                print("num_I_zero: ", len(self.i_zero_coords_list))
+                i_arr_zero = np.array(self.i_zero_coords_list) / ang_borr
+                i_hist_zero, _ = np.histogramdd(
+                    i_arr_zero, bins=self.mesh, range=bounds)
+                self.i_hist_zero = i_hist_zero.ravel() / total_I_num / volume
+
+            # I (Corner)
+            if len(self.i_corner_coords_list) > 0:
+                print("num_I_corner: ", len(self.i_corner_coords_list))
+                i_arr_corner = np.array(self.i_corner_coords_list) / ang_borr
+                i_hist_corner, _ = np.histogramdd(
+                    i_arr_corner, bins=self.mesh, range=bounds)
+                self.i_hist_corner = i_hist_corner.ravel() / total_I_num / volume
+            # I (Edge)
+            if len(self.i_edge_coords_list) > 0:
+                print("num_I_edge: ", len(self.i_edge_coords_list))
+                i_arr_edge = np.array(self.i_edge_coords_list) / ang_borr
+                i_hist_edge, _ = np.histogramdd(
+                    i_arr_edge, bins=self.mesh, range=bounds)
+                self.i_hist_edge = i_hist_edge.ravel() / total_I_num / volume
+                
+            # I (Three)
+            if len(self.i_three_coords_list) > 0:
+                print("num_I_three: ", len(self.i_three_coords_list))
+                i_arr_three = np.array(self.i_three_coords_list) / ang_borr
+                i_hist_three, _ = np.histogramdd(
+                    i_arr_three, bins=self.mesh, range=bounds)
+                self.i_hist_three = i_hist_three.ravel() / total_I_num / volume
+
+    def outputCube(self):
+        output_start_time = time.time()
+        base = re.match(r"(\d{3}.*?LiI).*?",
+                        args.trjfile.split("/")[-1]).group(1)
+        self.atomsDic = {'I': '53', 'Li': '3', 'P': '15', 'S': '16'}
+
+        def get_header():
+            body = f"created from {__file__}, {args}\n"
+            body += "Contains the selected quantity on a FFT grid\n"
+            origin = [-(self.cutoff)/ang_borr] * 3
+            body += "{:>5d}{:>12.7f}{:>12.7f}{:>12.7f}\n".format(5, *origin)
+            body += "{:>5d}{:>12.7f}{:>12.7f}{:>12.7f}\n".format(
+                self.mesh[0], (self.cutoff*2)/self.mesh[0]/ang_borr, 0, 0)
+            body += "{:>5d}{:>12.7f}{:>12.7f}{:>12.7f}\n".format(
+                self.mesh[1], 0, (self.cutoff*2)/self.mesh[1]/ang_borr, 0)
+            body += "{:>5d}{:>12.7f}{:>12.7f}{:>12.7f}\n".format(
+                self.mesh[2], 0, 0, (self.cutoff*2)/self.mesh[2]/ang_borr)
+            body += "{:>5d}{:>12.7f}{:>12.7f}{:>12.7f}{:>12.7f}\n".format(
+                int(self.atomsDic["P"]), float(self.atomsDic["P"]), *self.ps4_coord[0])
+            for s_coord in self.ps4_coord[1:5]:
+                body += "{:>5d}{:>12.7f}{:>12.7f}{:>12.7f}{:>12.7f}\n".format(
+                    int(self.atomsDic["S"]), float(self.atomsDic["S"]), *s_coord)
+            return body
+
+        def save_cube(hist_data, suffix):
+            if hist_data is None:
+                return
+            body = get_header()
+            for idx, r in enumerate(hist_data):
+                if idx % 6 == 5:
+                    body += "{:>13.5E}\n".format(r)
+                else:
+                    body += "{:>13.5E}".format(r)
+            if idx % 6 != 5:
+                body += "\n"  # 最後の行で改行がない場合用
+
+            outfile = f"{dirname}/{base}_PS4_{suffix}.cube"
+            with open(outfile, "w") as o:
+                o.write(body)
+            print(f"{outfile} was created.")
+
+        save_cube(self.li_hist, "Li")
+        save_cube(self.i_hist, "I_All")
+        save_cube(self.i_hist_zero, "I_zero")
+        save_cube(self.i_hist_corner, "I_corner")
+        save_cube(self.i_hist_edge, "I_edge")
+        save_cube(self.i_hist_three, "I_three")
+
+        output_end_time = time.time()
+        print(f"output_time : {output_end_time - output_start_time:.2f} s")
+
+
+if __name__ == "__main__":
+    trj = LammpsTrj()
+    trj.makeCube(args.trjfile)
+    trj.outputCube()