diff --git a/README.md b/README.md
index 689ab89..7e5d9e1 100644
--- a/README.md
+++ b/README.md
@@ -90,20 +90,22 @@ There are two options, <u>*interactive mode*</u> and <u>*command-line mode*</u>
 3. Follow the on-screen prompts to interactively select and run the desired function.
 
     ```
-             ____ ____  _   _ __  __ ____  _    _ _
-            / ___|  _ \| | | |  \/  |  _ \| | _(_) |_
-           | |  _| |_) | | | | |\/| | | | | |/ / | __|
-           | |_| |  __/| |_| | |  | | |_| |   <| | |_
-            \____|_|    \___/|_|  |_|____/|_|\_\_|\__|
+               ____ ____  _   _ __  __ ____  _    _ _
+              / ___|  _ \| | | |  \/  |  _ \| | _(_) |_
+             | |  _| |_) | | | | |\/| | | | | |/ / | __|
+             | |_| |  __/| |_| | |  | | |_| |   <| | |_
+              \____|_|    \___/|_|  |_|____/|_|\_\_|\__|
     
-            GPUMDkit Version 1.5.3 (dev) (2026-03-18)
-      Core Developer: Zihan YAN (yanzihan@westlake.edu.cn)
+              GPUMDkit Version 1.5.5 (dev) (2026-05-10)
+        Core Developer: Zihan YAN (yanzihan@westlake.edu.cn)
+     Main Contributors: Denan LI, Xin WU, Zhoulin LIU & Chen HUA
     
-     ----------------------- GPUMD -----------------------
+     ---------------------- GPUMD ------------------------
      1) Format Conversion          2) Sample Structures
      3) Workflow                   4) Calculators
-     5) Analyzer                   6) Developing ...
-     0) Quit!
+     5) Analyzer                   6) Visualization
+     7) Utilities                  8) Developing...
+     0) Exit
      ------------>>
      Input the function number:
     ```
@@ -123,29 +125,38 @@ gpumdkit.sh -h
 the help information:
 
 ```
-+==================================================================================================+
-|                              GPUMDkit 1.5.3 (dev) (2026-03-18) Usage                             |
-+======================================== Conversions =============================================+
-| -out2xyz       Convert OUTCAR to extxyz       | -pos2exyz     Convert POSCAR to extxyz           |
-| -cif2pos       Convert cif to POSCAR          | -pos2lmp      Convert POSCAR to LAMMPS           |
-| -cif2exyz      Convert cif to extxyz          | -lmp2exyz     Convert LAMMPS-dump to extxyz      |
-| -addgroup      Add group label                | -addweight    Add weight to the struct in extxyz |
-| -cp2k2xyz      Convert CP2K file to extxyz    | -traj2exyz    Convert ASE traj to extxyz         |
-| -xdat2exyz     Convert XDATCAR to extxyz      | Developing...                                    |
-+========================================= Analysis ===============================================+
-| -range         Print range of energy etc.     | -max_rmse     Get max RMSE from extxyz           |
-| -min_dist      Get min_dist between atoms     | -min_dist_pbc Get min_dist considering PBC       |
-| -filter_box    Filter struct by box limits    | -filter_value Filter struct by value (efs)       |
-| -filter_dist   Filter struct by min_dist      | -analyze_comp Analyze composition of extxyz      |
-| -pynep         Sample struct by pynep         | Developing...                                    |
-+====================================== Misc Utilities ============================================+
-| -plt           Plot scripts                   | -get_frame     Extract the specified frame       |
-| -calc          Calculators                    | -frame_range   Extract frames by fraction range  |
-| -clean         Clear files for work_dir       | -clean_xyz     Clean extra info in XYZ file      |
-| -time          Time consuming Analyzer        | -update        Update GPUMDkit                   |
-+==================================================================================================+
-| For detailed usage and examples, use: gpumdkit.sh -<option> -h                                   |
-+==================================================================================================+
++-------------------------------------------------------------------------------------------------------+
+|                          GPUMDkit 1.5.5 (dev) (2026-05-10) Command Help                               |
++-------------------------------------------------------------------------------------------------------+
+|                                          MAIN FUNCTIONS                                               |
++-------------------------------------------------------------------------------------------------------+
+| -h            Show this help table            | -plt <type>        Plot and visualization tools       |
+| -calc <type>  Calculator tools                | -time <gpumd|nep>  Time-consuming analyzer            |
+| -update       Update GPUMDkit                 | -clean             Clean extra files in current dir   |
++-------------------------------------------------------------------------------------------------------+
+|                                         FORMAT CONVERSION                                             |
++-------------------------------------------------------------------------------------------------------+
+| -out2xyz      OUTCAR -> extxyz (shell)        | -out2exyz          OUTCAR -> extxyz (python)          |
+| -cp2k2xyz     CP2K log -> xyz                 | -xdat2exyz         XDATCAR -> extxyz                  |
+| -cif2pos      cif -> POSCAR                   | -cif2exyz          cif -> extxyz                      |
+| -pos2exyz     POSCAR -> extxyz                | -exyz2pos          extxyz -> POSCAR                   |
+| -pos2lmp      POSCAR -> LAMMPS data           | -lmp2exyz          LAMMPS dump -> extxyz              |
+| -traj2exyz    ASE traj -> extxyz              | -replicate         Replicate structure                |
+| -addgroup     Add group labels                | -addweight         Add structure weight in extxyz     |
+| -clean_xyz    Clean extra info in extxyz      | -get_frame         Extract specific frame             |
+| -frame_range  Extract frames by range         |                                                       |
++-------------------------------------------------------------------------------------------------------+
+|                                            ANALYSIS                                                   |
++-------------------------------------------------------------------------------------------------------+
+| -range        Energy/force/virial statistics  | -analyze_comp      Analyze composition                |
+| -chem_species Analyze chemical species        | -cbc               Charge balance check               |
+| -min_dist     Min distance (no PBC)           | -min_dist_pbc      Min distance with PBC              |
+| -filter_dist  Filter by min_dist (no PBC)     | -filter_dist_pbc   Filter by min_dist (PBC)           |
+| -pda          Probability density analysis    | -hbond             Hydrogen-bond analysis             |
+| -pynep        FPS sampling by PyNEP           |                                                       |
++-------------------------------------------------------------------------------------------------------+
+| Detailed usage: gpumdkit.sh -<option> -h    Plot details: gpumdkit.sh -plt <type> -h                  |
++-------------------------------------------------------------------------------------------------------+
 ```
 
 ##### Example 2: View help information for -plt
@@ -157,28 +168,45 @@ gpumdkit.sh -plt -h
 the help information:
 
 ```
-+=====================================================================================================+
-|                              GPUMDkit 1.5.3 (dev) (2026-03-18) Plotting Usage                       |
-+=============================================== Plot Types ==========================================+
-| thermo          Plot thermo info                   | train          Plot NEP train results          |
-| prediction      Plot NEP prediction results        | train_test     Plot NEP train and test results |
-| msd             Plot mean square displacement      | msd_conv       Plot the convergence of MSD     |
-| msd_all         Plot MSD of all species            | sdc            Plot self diffusion coefficient |
-| rdf             Plot radial distribution function  | vac            Plot velocity autocorrelation   |
-| restart         Plot parameters in nep.restart     | dimer          Plot dimer plot                 |
-| force_errors    Plot force errors                  | des            Plot descriptors                |
-| charge          Plot charge distribution           | lr             Plot learning rate              |
-| doas            Plot density of atomistic states   | net_force      Plot net force distribution     |
-| sigma           Plot Arrhenius sigma               | D              Plot Arrhenius diffusivity      |
-| sigma_xyz       Plot directional Arrhenius sigma   | D_xyz          Plot directional Arrhenius D    |
-| emd             Plot EMD results                   | nemd           Plot NEMD results               |
-| hnemd           Plot HNEMD results                 | pdos           Plot VAC and PDOS               |
-| plane-grid      Plot displacement plane grid       | parity_density Plot parity plot density        |
-| cohesive        Plot cohsive energy                | viscosity      Plot visconsity                 |
-| rdf_pmf         Plot potential of mean force (PMF) |                                                |
-+=====================================================================================================+
-| For detailed usage and examples, use: gpumdkit.sh -plt <plot_type> -h                               |
-+=====================================================================================================+
+ +-----------------------------------------------------------------------------------------------+
+ |                     GPUMDkit 1.5.5 (dev) (2026-05-10) PLOT & VISUALIZATION TOOLS              |
+ +-----------------------------------------------------------------------------------------------+
+ |  Usage: gpumdkit.sh -plt <type>                        Help: gpumdkit.sh -plt <type> -h       |
+ +-----------------------------------------------------------------------------------------------+
+ |                                    NEP Training & Evaluation                                  |
+ +-----------------------------------------------------------------------------------------------+
+ |  train          - NEP training results           prediction     - NEP prediction results      |
+ |  train_test     - NEP train and test results     parity_density - Parity density plot         |
+ |  train_density  - Training results density plot  restart        - Parameters in nep.restart   |
+ |  charge         - Charge distribution            born_charge    - Born effective charges      |
+ |  dimer          - Dimer energy/force curve       force_errors   - Force errors                |
+ |  des            - Descriptors                    lr             - Learning rate for gnep      |
+ +-----------------------------------------------------------------------------------------------+
+ |                                     Diffusion & Transport                                     |
+ +-----------------------------------------------------------------------------------------------+
+ |  msd            - Mean square displacement       msd_conv       - MSD convergence             |
+ |  msd_all        - MSD for all species            sdc            - Self diffusion coefficient  |
+ |  msd_sdc        - MSD and SDC together           sigma          - Arrhenius ionic conductivity|
+ |  D              - Arrhenius diffusivity          sigma_xyz      - Directional Arrhenius sigma |
+ |  D_xyz          - Directional Arrhenius D                                                     |
+ +-----------------------------------------------------------------------------------------------+
+ |                                    MD & Structural Analysis                                   |
+ +-----------------------------------------------------------------------------------------------+
+ |  thermo         - thermo info in thermo.out      thermo2/3      - Thermo in different styles  |
+ |  rdf            - Radial distribution function   rdf_pmf        - Potential of mean force     |
+ |  vac            - Velocity autocorrelation       cohesive       - Cohesive energy curve       |
+ |  net_force      - Net force distribution         plane-grid     - Displacement plane grid     |
+ |  doas           - Density of atomistic states                                                 |
+ +-----------------------------------------------------------------------------------------------+
+ |                                        Heat Transport                                         |
+ +-----------------------------------------------------------------------------------------------+
+ |  emd            - EMD results                    nemd           - NEMD results                |
+ |  hnemd          - HNEMD results                  viscosity      - Viscosity                   |
+ +-----------------------------------------------------------------------------------------------+
+ |                                          Phonons                                              |
+ +-----------------------------------------------------------------------------------------------+
+ |  pdos           - VAC and PDOS                                                                |
+ +-----------------------------------------------------------------------------------------------+
 ```
 
 ##### Example 3: Convert VASP OUTCARs to extxyz
@@ -203,7 +231,19 @@ gpumdkit.sh -plt train
     <img src="./docs/Gallery/train.png" alt="msd" width="75%" />
 </div>
 
-##### Example 5: Plot thermo evolution
+##### Example 5: Plot the parity plots
+
+To visualize the parity plots:
+
+```
+gpumdkit.sh -plt test
+```
+
+<div align="center">
+    <img src="./docs/Gallery/prediction.png" alt="msd" width="95%" />
+</div>
+
+##### Example 6: Plot thermo evolution
 
 To visualize `thermo` evolution from `thermo.out` :
 
diff --git a/Scripts/calculators/calc_minimize.py b/Scripts/calculators/calc_minimize.py
new file mode 100644
index 0000000..d48e18d
--- /dev/null
+++ b/Scripts/calculators/calc_minimize.py
@@ -0,0 +1,89 @@
+"""
+This script performs structure relaxation using calorine.
+It reads a single structure from POSCAR or extxyz file, performs energy minimization,
+and outputs the optimization trajectory to minimize.xyz in extxyz format.
+The force convergence threshold and maximum steps can be specified via command line arguments.
+"""
+
+import sys
+from tqdm import tqdm
+import numpy as np
+from ase import Atoms
+from ase.io import read, write
+from ase.optimize import BFGS
+from calorine.calculators import CPUNEP
+import os
+
+class TrajectoryCallback:
+    """Callback class to save trajectory during optimization."""
+    def __init__(self, atoms):
+        self.atoms = atoms
+        self.trajectory = [atoms.copy()]
+    
+    def __call__(self):
+        self.trajectory.append(self.atoms.copy())
+
+def main():
+    # Parse command line arguments
+    if len(sys.argv) < 3:
+        print(" Usage: python calc_minimize.py <structure_file> <nep.txt> [fmax=0.01] [max_steps=1000]")
+        print(" Supported structure formats: POSCAR, .xyz")
+        print(" Optional arguments: ")
+        print("  fmax: Force convergence threshold in eV/Ang (default: 0.01)")
+        print("  max_steps: Maximum number of optimization steps (default: 1000)")
+        sys.exit(1)
+
+    structure_file = sys.argv[1]  # Structure file (POSCAR or extxyz)
+    model_path = sys.argv[2]  # Path to the model for CPUNEP
+
+    # Check if structure file exists
+    if not os.path.exists(structure_file):
+        print(f" Error: Structure file '{structure_file}' does not exist.")
+        sys.exit(1)
+
+    # Check if model file exists
+    if not os.path.exists(model_path):
+        print(f" Error: Model file '{model_path}' does not exist.")
+        sys.exit(1)
+
+    # Set default values and parse optional arguments
+    fmax = float(sys.argv[3]) if len(sys.argv) > 3 else 0.01  # Default max force threshold
+    max_steps = int(sys.argv[4]) if len(sys.argv) > 4 else 1000  # Default max steps
+
+    print(f" Force convergence threshold: {fmax} eV/Ang")
+    print(f" Maximum steps: {max_steps}")
+    print(f" Structure file: {structure_file}")
+    print(f" Model path: {model_path}")
+
+    # Determine file format and read structure
+    file_ext = os.path.splitext(structure_file)[1].lower()
+    if file_ext in ['.xyz']:
+        atoms = read(structure_file)
+    elif 'POSCAR' in structure_file or 'CONTCAR' in structure_file:
+        atoms = read(structure_file, format='vasp')
+    else:
+        print(f" Warning: Unknown file extension '{file_ext}'. Attempting to read as extxyz.")
+        atoms = read(structure_file)
+
+    print(f" Read structure with {len(atoms)} atoms")
+
+    # Create calculator and assign to atoms
+    calc = CPUNEP(model_path)
+    atoms.calc = calc
+
+    # Initialize trajectory callback
+    traj_callback = TrajectoryCallback(atoms)
+
+    # Perform structure optimization with trajectory tracking
+    opt = BFGS(atoms, logfile='minimize.log')
+    opt.attach(traj_callback)  # Attach callback to save trajectory
+    opt.run(fmax=fmax, steps=max_steps)  # Convergence criterion based on input
+
+    # Save the complete trajectory to minimize.xyz
+    write('minimize.xyz', traj_callback.trajectory)
+    print(f" Optimization completed after {len(traj_callback.trajectory)} steps.")
+    print(f" Final max force: {np.max(np.abs(atoms.get_forces())):.6f} eV/Ang")
+    print(f" Optimization trajectory saved to minimize.xyz")
+
+if __name__ == "__main__":
+    main()
diff --git a/Scripts/calculators/calc_msd.py b/Scripts/calculators/calc_msd.py
new file mode 100644
index 0000000..4eb0af0
--- /dev/null
+++ b/Scripts/calculators/calc_msd.py
@@ -0,0 +1,405 @@
+"""
+Mean Square Displacement (MSD) Calculator from extxyz Trajectory
+
+Reads an extxyz trajectory file, extracts positions of a target element,
+unwraps coordinates across periodic boundaries if needed, and computes
+the directional MSD (x, y, z) via the Wiener-Khinchin FFT algorithm.
+
+Usage:
+    python calc_msd.py <extxyz_file> <element_symbol> <dt_fs> [max_corr_steps]
+
+Arguments:
+    extxyz_file     Path to the input extxyz trajectory file
+    element_symbol  Chemical symbol of the element to analyze (e.g. Li, O, Na)
+    dt_fs           Time step between consecutive frames [fs]
+    max_corr_steps  (Optional) Maximum number of correlation lag steps to output.
+                    Output time range will be [dt_fs, max_corr_steps * dt_fs].
+                    If omitted, all lag times up to (n_frames - 2) are written.
+
+Output:
+    msd.out - Four-column ASCII file:
+              Col 1: Time (ps)
+              Col 2: MSD_x (Ang^2)
+              Col 3: MSD_y (Ang^2)
+              Col 4: MSD_z (Ang^2)
+
+Example:
+    python calc_msd.py traj.xyz Li 2.0
+    python calc_msd.py traj.xyz Li 2.0 5000
+"""
+
+import sys
+import os
+import time
+import numpy as np
+
+
+# ==============================================================================
+# Argument parsing
+# ==============================================================================
+
+def parse_args():
+    """
+    Parse and validate positional command-line arguments.
+
+    Returns
+    -------
+    extxyz_file    : str
+    element_symbol : str
+    dt_fs          : float
+    max_corr_steps : int or None
+    """
+    args = sys.argv[1:]
+
+    if len(args) == 0 or args[0] in ("-h", "--help"):
+        print(__doc__)
+        sys.exit(0)
+
+    if len(args) not in (3, 4):
+        print(" Error: Expected 3 or 4 positional arguments, got %d." % len(args))
+        print("        Run with -h for usage information.")
+        sys.exit(1)
+
+    extxyz_file    = args[0]
+    element_symbol = args[1]
+    dt_str         = args[2]
+    max_corr_str   = args[3] if len(args) == 4 else None
+
+    # Validate file
+    if not os.path.isfile(extxyz_file):
+        print(" Error: File not found: '%s'" % extxyz_file)
+        sys.exit(1)
+
+    if not extxyz_file.lower().endswith((".xyz", ".extxyz")):
+        print(" Warning: '%s' does not have a recognised .xyz / .extxyz extension. "
+              " Proceeding anyway." % extxyz_file)
+
+    # Validate element symbol
+    if not element_symbol.isalpha() or not element_symbol[0].isupper():
+        print(" Warning: '%s' may not be a valid element symbol. "
+              " Expected format: e.g. 'Li', 'O', 'Na'." % element_symbol)
+
+    # Validate dt
+    try:
+        dt_fs = float(dt_str)
+    except ValueError:
+        print(" Error: dt_fs must be a number, got: '%s'" % dt_str)
+        sys.exit(1)
+    if dt_fs <= 0:
+        print(" Error: dt_fs must be positive, got: %s" % dt_fs)
+        sys.exit(1)
+
+    # Validate max_corr_steps (optional)
+    max_corr_steps = None
+    if max_corr_str is not None:
+        try:
+            max_corr_steps = int(max_corr_str)
+        except ValueError:
+            print(" Error: max_corr_steps must be a positive integer, got: '%s'" % max_corr_str)
+            sys.exit(1)
+        if max_corr_steps <= 0:
+            print(" Error: max_corr_steps must be a positive integer, got: %d" % max_corr_steps)
+            sys.exit(1)
+
+    return extxyz_file, element_symbol, dt_fs, max_corr_steps
+
+
+# ==============================================================================
+# Trajectory I/O
+# ==============================================================================
+
+def read_extxyz_file(filename, element_symbol):
+    """
+    Read coordinates of the specified element from an extxyz file.
+
+    Parameters
+    ----------
+    filename       : str
+    element_symbol : str
+
+    Returns
+    -------
+    all_positions : ndarray, shape (n_frames, n_atoms, 3)
+    cell_info     : ndarray, shape (3, 3), or None
+    """
+    try:
+        from ase.io import read
+    except ImportError:
+        print(" Error: ASE is required but not installed.")
+        print("        Install it with: pip install ase")
+        sys.exit(1)
+
+    try:
+        atoms_list = read(filename, index=":")
+    except Exception as exc:
+        print(" Error: Failed to read '%s': %s" % (filename, exc))
+        sys.exit(1)
+
+    if not atoms_list:
+        print(" Error: No frames found in '%s'." % filename)
+        sys.exit(1)
+
+    all_positions = []
+    cell_info = None
+
+    for frame_idx, atoms in enumerate(atoms_list):
+        if frame_idx == 0:
+            cell_info = np.array(atoms.get_cell())
+
+        symbols = atoms.get_chemical_symbols()
+        element_indices = [i for i, s in enumerate(symbols) if s == element_symbol]
+
+        if not element_indices:
+            print(" Error: Element '%s' not found in frame %d. "
+                  " Check the element symbol and trajectory file."
+                  % (element_symbol, frame_idx))
+            sys.exit(1)
+
+        all_positions.append(atoms.positions[element_indices])
+
+    all_positions = np.array(all_positions)   # shape: (n_frames, n_atoms, 3)
+    return all_positions, cell_info
+
+
+# ==============================================================================
+# Coordinate wrapping utilities
+# ==============================================================================
+
+def check_wrapped_coordinates(positions, cell):
+    """
+    Heuristically check whether Cartesian coordinates are PBC-wrapped.
+
+    Returns True if any coordinate is negative, exceeds 1.5x the cell
+    length, or if the per-dimension coordinate span exceeds the cell length.
+    """
+    if cell is None or positions.size == 0:
+        return False
+
+    cell_lengths = np.array([np.linalg.norm(cell[i]) for i in range(3)])
+
+    for dim in range(3):
+        col = positions[:, :, dim]
+        if np.any(col < 0) or np.any(col > cell_lengths[dim] * 1.5):
+            return True
+        if col.max() - col.min() > cell_lengths[dim] * 1.5:
+            return True
+
+    return False
+
+
+def unwrap_coordinates(positions, cell):
+    """
+    Unwrap PBC-wrapped coordinates using the minimum-image convention.
+
+    For each atom and each consecutive frame pair, if the displacement
+    along any axis exceeds half the cell length, all subsequent frames
+    are shifted by one full cell vector to remove the jump.
+
+    Parameters
+    ----------
+    positions : ndarray, shape (n_frames, n_atoms, 3)
+    cell      : ndarray, shape (3, 3)
+
+    Returns
+    -------
+    ndarray, shape (n_frames, n_atoms, 3)
+    """
+    if positions.shape[0] < 2 or cell is None:
+        return positions
+
+    cell_lengths = np.array([np.linalg.norm(cell[i]) for i in range(3)])
+    unwrapped = np.copy(positions)
+
+    for atom_idx in range(positions.shape[1]):
+        for frame in range(1, positions.shape[0]):
+            displacement = positions[frame, atom_idx] - positions[frame - 1, atom_idx]
+            for dim in range(3):
+                half = cell_lengths[dim] / 2.0
+                if displacement[dim] > half:
+                    unwrapped[frame:, atom_idx, dim] -= cell_lengths[dim]
+                elif displacement[dim] < -half:
+                    unwrapped[frame:, atom_idx, dim] += cell_lengths[dim]
+
+    return unwrapped
+
+
+# ==============================================================================
+# MSD calculation
+# ==============================================================================
+
+def msd_fft(positions):
+    """
+    Compute directional MSD per axis via the Wiener-Khinchin FFT algorithm.
+
+    Uses the identity:
+        MSD(dt) = 2*<r^2> - 2*C(dt)
+    where C is the position autocorrelation function evaluated via FFT,
+    giving O(N log N) complexity.
+
+    Parameters
+    ----------
+    positions : ndarray, shape (n_frames, n_atoms, 3)
+
+    Returns
+    -------
+    msd_xyz : ndarray, shape (n_frames - 2, 3)
+        Columns correspond to MSD_x, MSD_y, MSD_z, each averaged over all
+        atoms. Lag index starts from 1 (first row = dt * 1 frame).
+    """
+    axis_time = 0
+    axis_atom = 1
+    n_time    = positions.shape[axis_time]
+    sq_pos    = np.square(positions)
+
+    # Autocorrelation via FFT (Wiener-Khinchin theorem)
+    fft_pos = np.fft.fft(positions, n=2 * n_time, axis=axis_time)
+    S2 = np.fft.ifft(np.abs(fft_pos) ** 2, axis=axis_time).real
+    S2 = S2[:n_time]                              # keep first half only
+
+    # Direct-sum (S1) term
+    sq_aug   = np.concatenate(
+        [sq_pos, np.zeros((1,) + sq_pos.shape[1:])], axis=axis_time
+    )                                             # shape: (n_time+1, n_atoms, 3)
+    total    = 2.0 * np.sum(sq_aug, axis=axis_time, keepdims=True)
+    forward  = np.concatenate(
+        [np.zeros((1,) + sq_pos.shape[1:]), sq_pos], axis=axis_time
+    )
+    backward = np.flip(sq_aug, axis=axis_time)
+    S1       = (total - np.cumsum(forward + backward, axis=axis_time))[:n_time]
+
+    # Per-atom MSD: shape (n_time, n_atoms, 3)
+    counts   = (n_time - np.arange(n_time)).reshape(-1, 1, 1)
+    atom_msd = (S1 - 2.0 * S2) / counts
+
+    # Drop lag=0 and the last lag (statistically unreliable)
+    lag_range = np.arange(1, n_time - 1)
+    atom_msd  = atom_msd[lag_range]              # shape: (n_time-2, n_atoms, 3)
+
+    # Average over atoms, keep x/y/z separate -> shape: (n_time-2, 3)
+    msd_xyz = atom_msd.mean(axis=axis_atom)
+    return msd_xyz
+
+
+# ==============================================================================
+# Main workflow
+# ==============================================================================
+
+def calculate_msd_from_extxyz(extxyz_file, element_symbol, dt_fs, max_corr_steps):
+    """
+    Full MSD calculation pipeline:
+      1. Read trajectory
+      2. Check / unwrap coordinates
+      3. Compute directional MSD via FFT
+      4. Truncate to max_corr_steps if specified
+      5. Write msd.out (time, MSD_x, MSD_y, MSD_z)
+    """
+    print(" ")
+    print(" "+ "=" * 60)
+    print("  Trajectory file  : %s" % extxyz_file)
+    print("  Target element   : %s" % element_symbol)
+    print("  Time step (dt)   : %.4f fs" % dt_fs)
+    if max_corr_steps is not None:
+        print("  Max corr. steps  : %d  (%.4f ps)"
+              % (max_corr_steps, max_corr_steps * dt_fs / 1000.0))
+    else:
+        print("  Max corr. steps  : all frames")
+    print(" "+ "-" * 60)
+
+    # ------------------------------------------------------------------
+    # Read trajectory
+    # ------------------------------------------------------------------
+    print(" ")
+    print(" Reading file: %s" % extxyz_file)
+    all_positions, cell = read_extxyz_file(extxyz_file, element_symbol)
+    n_frames, n_atoms, _ = all_positions.shape
+    print("   Loaded %d frames, %d %s atoms per frame"
+          % (n_frames, n_atoms, element_symbol))
+
+    if cell is not None:
+        cell_lengths = np.array([np.linalg.norm(cell[i]) for i in range(3)])
+        print("   Cell lengths (a, b, c) : %.4f  %.4f  %.4f  Ang"
+              % (cell_lengths[0], cell_lengths[1], cell_lengths[2]))
+    else:
+        print(" Warning: No cell information found. PBC unwrapping will be skipped.")
+
+    # ------------------------------------------------------------------
+    # Coordinate wrapping check
+    # ------------------------------------------------------------------
+    print(" ")
+    # print(" "+ "-" * 60)
+    if cell is not None:
+        is_wrapped = check_wrapped_coordinates(all_positions, cell)
+        if is_wrapped:
+            print(" Coordinate status : Wrapped (PBC detected)")
+            print(" Unwrapping coordinates...")
+            all_positions = unwrap_coordinates(all_positions, cell)
+            print(" Coordinate unwrapping completed.")
+        else:
+            print(" Coordinate status : Unwrapped (no action needed)")
+    else:
+        print(" Coordinate status : Skipped (no cell information)")
+
+    # ------------------------------------------------------------------
+    # MSD calculation
+    # ------------------------------------------------------------------
+    print(" ")
+    # print(" "+ "-" * 60)
+    if n_frames < 3:
+        print(" Error: At least 3 frames are required for MSD, "
+              " but only %d found." % n_frames)
+        sys.exit(1)
+
+    print(" Calculating MSD (Wiener-Khinchin FFT)...")
+    msd_xyz = msd_fft(all_positions)             # shape: (n_frames-2, 3)
+
+    # Build lag-time axis: lag 1, 2, ..., n_frames-2
+    n_lags  = msd_xyz.shape[0]
+    time_ps = np.arange(1, n_lags + 1) * dt_fs / 1000.0   # fs -> ps
+
+    # Truncate to max_corr_steps if requested
+    if max_corr_steps is not None:
+        if max_corr_steps > n_lags:
+            print(" Warning: max_corr_steps (%d) exceeds available lag steps (%d). "
+                  " Using all available steps." % (max_corr_steps, n_lags))
+        else:
+            time_ps = time_ps[:max_corr_steps]
+            msd_xyz = msd_xyz[:max_corr_steps]
+
+    print("   MSD computed for %d lag times" % len(time_ps))
+    print("   Lag-time range  : %.6f ps  ->  %.6f ps"
+          % (time_ps[0], time_ps[-1]))
+
+    # ------------------------------------------------------------------
+    # Write output
+    # ------------------------------------------------------------------
+    print(" ")
+    # print(" "+ "-" * 60)
+    output_file = "msd.out"
+    header = ("msd.out by GPUMDkit\n"
+               "Element: %s  |  dt: %.4f fs  |  Frames: %d  |  Atoms: %d\n"
+               "%12s  %14s  %14s  %14s"
+               % (element_symbol, dt_fs, n_frames, n_atoms,
+                  "Time(ps)", "MSD_x(Ang^2)", "MSD_y(Ang^2)", "MSD_z(Ang^2)"))
+
+    np.savetxt(
+        output_file,
+        np.column_stack((time_ps, msd_xyz)),
+        fmt="%14.6f",
+        header=header,
+    )
+
+    print(" MSD calculation completed. Results saved to '%s'." % output_file)
+    print(" You can use [gpumdkit.sh -plt msd] to visualize the results.")
+
+# ==============================================================================
+# Entry point
+# ==============================================================================
+
+if __name__ == "__main__":
+    t0 = time.perf_counter()
+
+    extxyz_file, element_symbol, dt_fs, max_corr_steps = parse_args()
+    calculate_msd_from_extxyz(extxyz_file, element_symbol, dt_fs, max_corr_steps)
+
+    print(" Total elapsed time: %.2f s" % (time.perf_counter() - t0))
+    print(" "+ "=" * 60)
\ No newline at end of file
diff --git a/Scripts/plt_scripts/README.md b/Scripts/plt_scripts/README.md
index ad47aec..7d50d14 100644
--- a/Scripts/plt_scripts/README.md
+++ b/Scripts/plt_scripts/README.md
@@ -252,8 +252,8 @@ Plots MSD for all atomic species separately when using `all_groups` in GPUMD.
 
 **Usage:**
 ```bash
-gpumdkit.sh -plt msd_all <species1> <species2> ...
-gpumdkit.sh -plt msd_all Li P S
+gpumdkit.sh -plt msd_all msd.out <species1> <species2> ...
+gpumdkit.sh -plt msd_all msd.out Li P S
 ```
 
 **Requirements:** Must use `all_groups` in the `compute_msd` command in `run.in`.
@@ -268,7 +268,7 @@ gpumdkit.sh -plt msd_all Li P S
 
 Checks convergence of MSD calculations across different time windows.
 
-**Input File:** `msd.out` (computed with `save_every` option)
+**Input File:** `msd_step*.out` (computed with `save_every` option)
 
 **Usage:**
 ```bash
@@ -663,7 +663,7 @@ gpumdkit.sh -plt net_force train.xyz
 | `nemd` | NEMD outputs | NEMD thermal transport |
 | `hnemd` | HNEMD outputs | HNEMD thermal transport |
 | `arrhenius_d` | Multiple MSD | Arrhenius diffusivity |
-| `arrhenius_sigma` | Conductivity data | Arrhenius conductivity |
+| `arrhenius_sigma` | Multiple MSD | Arrhenius conductivity |
 
 ## Contributing
 
@@ -681,26 +681,46 @@ See [CONTRIBUTING.md](../../CONTRIBUTING.md) for detailed guidelines.
 
 For quick reference, here's the complete command list:
 
-```bash
-gpumdkit.sh -plt
-+=====================================================================================================+
-|                              GPUMDkit 1.4.2 (dev) (2025-12-17) Plotting Usage                      |
-+=============================================== Plot Types ==========================================+
-| thermo          Plot thermo info                   | train          Plot NEP train results          |
-| prediction      Plot NEP prediction results        | train_test     Plot NEP train and test results |
-| msd             Plot mean square displacement      | msd_conv       Plot the convergence of MSD     |
-| msd_all         Plot MSD of all species            | sdc            Plot self diffusion coefficient |
-| rdf             Plot radial distribution function  | vac            Plot velocity autocorrelation   |
-| restart         Plot parameters in nep.restart     | dimer          Plot dimer plot                 |
-| force_errors    Plot force errors                  | des            Plot descriptors                |
-| charge          Plot charge distribution           | lr             Plot learning rate              |
-| doas            Plot density of atomistic states   | arrhenius_d    Plot Arrhenius diffusivity      |
-| arrhenius_sigma Plot Arrhenius sigma               | net_force      Plot net force distribution     |
-| emd             Plot EMD results                   | nemd           Plot NEMD results               |
-| hnemd           Plot HNEMD results                 | parity_density Parity plots with density       |
-+=====================================================================================================+
-| For detailed usage and examples, use: gpumdkit.sh -plt <plot_type> -h                              |
-+=====================================================================================================+
+```
+ +-----------------------------------------------------------------------------------------------+
+ |                     GPUMDkit 1.5.5 (dev) (2026-05-10) PLOT & VISUALIZATION TOOLS              |
+ +-----------------------------------------------------------------------------------------------+
+ |  Usage: gpumdkit.sh -plt <type>                        Help: gpumdkit.sh -plt <type> -h       |
+ +-----------------------------------------------------------------------------------------------+
+ |                                    NEP Training & Evaluation                                  |
+ +-----------------------------------------------------------------------------------------------+
+ |  train          - NEP training results           prediction     - NEP prediction results      |
+ |  train_test     - NEP train and test results     parity_density - Parity density plot         |
+ |  train_density  - Training results density plot  restart        - Parameters in nep.restart   |
+ |  charge         - Charge distribution            born_charge    - Born effective charges      |
+ |  dimer          - Dimer energy/force curve       force_errors   - Force errors                |
+ |  des            - Descriptors                    lr             - Learning rate for gnep      |
+ +-----------------------------------------------------------------------------------------------+
+ |                                     Diffusion & Transport                                     |
+ +-----------------------------------------------------------------------------------------------+
+ |  msd            - Mean square displacement       msd_conv       - MSD convergence             |
+ |  msd_all        - MSD for all species            sdc            - Self diffusion coefficient  |
+ |  msd_sdc        - MSD and SDC together           sigma          - Arrhenius ionic conductivity|
+ |  D              - Arrhenius diffusivity          sigma_xyz      - Directional Arrhenius sigma |
+ |  D_xyz          - Directional Arrhenius D                                                     |
+ +-----------------------------------------------------------------------------------------------+
+ |                                    MD & Structural Analysis                                   |
+ +-----------------------------------------------------------------------------------------------+
+ |  thermo         - thermo info in thermo.out      thermo2/3      - Thermo in different styles  |
+ |  rdf            - Radial distribution function   rdf_pmf        - Potential of mean force     |
+ |  vac            - Velocity autocorrelation       cohesive       - Cohesive energy curve       |
+ |  net_force      - Net force distribution         plane-grid     - Displacement plane grid     |
+ |  doas           - Density of atomistic states                                                 |
+ +-----------------------------------------------------------------------------------------------+
+ |                                        Heat Transport                                         |
+ +-----------------------------------------------------------------------------------------------+
+ |  emd            - EMD results                    nemd           - NEMD results                |
+ |  hnemd          - HNEMD results                  viscosity      - Viscosity                   |
+ +-----------------------------------------------------------------------------------------------+
+ |                                          Phonons                                              |
+ +-----------------------------------------------------------------------------------------------+
+ |  pdos           - VAC and PDOS                                                                |
+ +-----------------------------------------------------------------------------------------------+
 ```
 
 ---
diff --git a/Scripts/plt_scripts/plt_arrhenius_d.py b/Scripts/plt_scripts/plt_arrhenius_d.py
index 2ce5c34..c106978 100644
--- a/Scripts/plt_scripts/plt_arrhenius_d.py
+++ b/Scripts/plt_scripts/plt_arrhenius_d.py
@@ -1,12 +1,26 @@
 """
-This script reads msd.out files from temperature-specific folders (e.g., 600K, 700K, etc.), 
-calculates the diffusion coefficients, and plots them in an Arrhenius plot. 
-It also fits the data to extract the activation energy. 
-
-Author: Modified from original by Zihan YAN (yanzihan@westlake.edu.cn)
-
-Usage:
-    python plt_arrhenius_d.py [save]
+=============================================================================
+GPUMDkit: A User-Friendly Toolkit for GPUMD and NEP
+Repository: https://github.com/zhyan0603/GPUMDkit
+=============================================================================
+Script:     plt_arrhenius_d.py
+Category:   Plot Scripts
+Purpose:    Compute diffusion coefficients from msd.out in temperature folders
+            (e.g., 600K, 700K, etc.), generate an Arrhenius plot, and extract 
+            activation energy.
+
+Usage:      gpumdkit.sh -plt D/arrhenius_d [save]
+            python plt_arrhenius_d.py [save]
+
+Arguments:
+  save    Save the plot as 'Arrhenius_D.png' instead of displaying it
+
+Output:
+  Arrhenius_D.png  (if save is used, or if backend is non-interactive)
+
+Author:     Zihan YAN (yanzihan@westlake.edu.cn)
+Last-modified: 2026-05-15
+=============================================================================    
 """
 
 import os
diff --git a/Scripts/plt_scripts/plt_prediction.py b/Scripts/plt_scripts/plt_prediction.py
index aa27ae2..12dd13d 100644
--- a/Scripts/plt_scripts/plt_prediction.py
+++ b/Scripts/plt_scripts/plt_prediction.py
@@ -1,104 +1,260 @@
 import sys
-import matplotlib.pyplot as plt
 import numpy as np
-# from cycler import cycler
+import matplotlib
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+# If running on a server without GUI, keep this commented unless needed
+# matplotlib.use('Agg')
 
-# custom_colors = ['#8FA2CD', '#F8BC7E', '#A9CA70', '#F09BA0', '#9D9EA3', '#B7B7EB']
-# plt.rcParams['axes.prop_cycle'] = cycler(color=custom_colors)
+# =========================
+# Global style
+# =========================
+plt.rcParams['font.family'] = 'sans-serif'
+plt.rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans', 'Liberation Sans']
+plt.rcParams['axes.unicode_minus'] = False
 
-# Load data
+plt.rcParams['axes.labelsize'] = 13
+plt.rcParams['xtick.labelsize'] = 12
+plt.rcParams['ytick.labelsize'] = 12
+
+plt.rcParams['axes.linewidth'] = 1.5
+plt.rcParams['xtick.major.width'] = 1.2
+plt.rcParams['ytick.major.width'] = 1.2
+plt.rcParams['xtick.color'] = 'black'
+plt.rcParams['ytick.color'] = 'black'
+plt.rcParams['axes.labelcolor'] = 'black'
+plt.rcParams['text.color'] = 'black'
+
+# =========================
+# Load data - No longer loading loss.out
+# =========================
 energy_data = np.loadtxt('energy_train.out')
 force_data = np.loadtxt('force_train.out')
-#virial_data = np.loadtxt('virial_train.out')
 stress_data = np.loadtxt('stress_train.out')
 
-# Filter out rows with invalid stress data
-valid_rows = ~np.any(np.abs(stress_data[:, :12]) >= 1e6, axis=1)
+# Filter invalid stress rows
+valid_rows = ~np.any(np.abs(stress_data[:, :12]) > 1e6, axis=1)
 stress_data = stress_data[valid_rows]
 
-# Function to calculate RMSE
-def calculate_rmse(pred, actual):
-    return np.sqrt(np.mean((pred - actual) ** 2))
+# =========================
+# Utility functions
+# =========================
+def rmse(pred, true):
+    return np.sqrt(np.mean((pred - true) ** 2))
+
+def mae(pred, true):
+    return np.mean(np.abs(pred - true))
 
-# Function to calculate dynamic axis limits
-def calculate_limits(train_data, padding=0.08):
-    data_min = np.min(train_data)
-    data_max = np.max(train_data)
+def r2_score_np(true, pred):
+    true = np.asarray(true).reshape(-1)
+    pred = np.asarray(pred).reshape(-1)
+    ss_res = np.sum((true - pred) ** 2)
+    ss_tot = np.sum((true - np.mean(true)) ** 2)
+    if ss_tot == 0:
+        return np.nan
+    return 1.0 - ss_res / ss_tot
+
+def get_limits(true, pred, padding=0.05):
+    true = np.asarray(true).reshape(-1)
+    pred = np.asarray(pred).reshape(-1)
+    data_min = min(np.min(true), np.min(pred))
+    data_max = max(np.max(true), np.max(pred))
     data_range = data_max - data_min
-    return data_min - padding * data_range, data_max + padding * data_range
-
-# Create a subplot with 1 row and 3 columns
-fig, axs = plt.subplots(1, 3, figsize=(12, 3.3), dpi=100)
-
-# Plotting the train_data figure
-xmin_energy, xmax_energy = calculate_limits(energy_data[:, 1])
-axs[0].set_xlim(xmin_energy, xmax_energy)
-axs[0].set_ylim(xmin_energy, xmax_energy)
-axs[0].plot(energy_data[:, 1], energy_data[:, 0], '.', markersize=10)
-axs[0].plot([xmin_energy, xmax_energy], [xmin_energy, xmax_energy], linewidth=2, color='grey', linestyle='--')
-axs[0].set_xlabel('DFT energy (eV/atom)', fontsize=11)
-axs[0].set_ylabel('NEP energy (eV/atom)', fontsize=11)
-axs[0].legend(['energy'], loc='upper left')
-axs[0].axis('tight')
-axs[0].tick_params(axis='both', labelsize=11)
-
-# Calculate and display RMSE for energy
-energy_rmse = calculate_rmse(energy_data[:, 0], energy_data[:, 1]) * 1000
-axs[0].text(0.35, 0.08, f'RMSE: {energy_rmse:.2f} meV/atom', transform=axs[0].transAxes, fontsize=11, verticalalignment='center')
-#axs[0].text(-0.1, 1.03, "(a)", transform=axs[0].transAxes, fontsize=13, va='top', ha='right')
-
-# Plotting the force_data figure
-xmin_force, xmax_force = calculate_limits(force_data[:, 3:6].reshape(-1))
-axs[1].set_xlim(xmin_force, xmax_force)
-axs[1].set_ylim(xmin_force, xmax_force)
-axs[1].plot(force_data[:, 3:6], force_data[:, 0:3], '.', markersize=10)
-axs[1].plot([xmin_force, xmax_force], [xmin_force, xmax_force], linewidth=2, color='grey', linestyle='--')
-axs[1].set_xlabel(r'DFT force (eV/$\mathrm{\AA}$)', fontsize=11)
-axs[1].set_ylabel(r'NEP force (eV/$\mathrm{\AA}$)', fontsize=11)
-axs[1].tick_params(axis='both', labelsize=11)
-axs[1].legend(['fx', 'fy', 'fz'], loc='upper left')
-axs[1].axis('tight')
-
-# Calculate and display RMSE for forces
-force_pred_flat = force_data[:, 0:3].flatten()
-force_actual_flat = force_data[:, 3:6].flatten()
-force_rmse = calculate_rmse(force_pred_flat, force_actual_flat) * 1000
-mean_force_rmse = force_rmse
-axs[1].text(0.35, 0.08, rf'RMSE: {mean_force_rmse:.2f} meV/$\mathrm{{\AA}}$', transform=axs[1].transAxes, fontsize=11, verticalalignment='center')
-#axs[1].text(-0.1, 1.03, "(b)", transform=axs[1].transAxes, fontsize=13, va='top', ha='right')
-
-# Plotting the stress figure
-xmin_stress, xmax_stress = calculate_limits(stress_data[:, 6:12].reshape(-1))
-axs[2].set_xlim(xmin_stress, xmax_stress)
-axs[2].set_ylim(xmin_stress, xmax_stress)
-axs[2].plot(stress_data[:, 6:12], stress_data[:, 0:6], '.', markersize=10)
-axs[2].plot([xmin_stress, xmax_stress], [xmin_stress, xmax_stress], linewidth=2, color='grey', linestyle='--')
-axs[2].set_xlabel('DFT stress (GPa)', fontsize=11)
-axs[2].set_ylabel('NEP stress (GPa)', fontsize=11)
-axs[2].tick_params(axis='both', labelsize=11)
-axs[2].legend(['xx', 'yy', 'zz', 'xy', 'yz', 'zx'], loc='upper left')
-axs[2].axis('tight')
-
-# Calculate and display RMSE for stresses
-stress_pred_flat = stress_data[:, 0:6].flatten()
-stress_actual_flat = stress_data[:, 6:12].flatten()
-stress_rmse = calculate_rmse(stress_pred_flat, stress_actual_flat)
-mean_stress_rmse = stress_rmse
-axs[2].text(0.35, 0.08, f'RMSE: {mean_stress_rmse:.4f} GPa', transform=axs[2].transAxes, fontsize=11, verticalalignment='center')
-#axs[2].text(-0.1, 1.03, "(c)", transform=axs[2].transAxes, fontsize=13, va='top', ha='right')
-
-# Adjust layout for better spacing
+    if data_range == 0:
+        data_range = 1.0
+    pad = padding * data_range
+    return data_min - pad, data_max + pad
+
+def beautify_axes(ax):
+    # ax.grid(True, linestyle='-', linewidth=0.8, alpha=0.5, color='lightgray')
+    ax.spines['top'].set_visible(False)
+    ax.spines['right'].set_visible(False)
+    ax.tick_params(axis='both', colors='black')
+
+def add_marginal_distributions(ax, true, pred, color, bins=35):
+    divider = make_axes_locatable(ax)
+
+    ax_top = divider.append_axes("top", size="18%", pad=0.06, sharex=ax)
+    ax_right = divider.append_axes("right", size="18%", pad=0.06, sharey=ax)
+
+    ax_top.hist(true, bins=bins, color=color, alpha=0.55, edgecolor='gray', linewidth=0.6)
+    ax_right.hist(pred, bins=bins, orientation='horizontal',
+                  color=color, alpha=0.55, edgecolor='gray', linewidth=0.6)
+
+    # Top marginal
+    ax_top.tick_params(axis='x', labelbottom=False, bottom=False)
+    ax_top.tick_params(axis='y', left=False, labelleft=False)
+    ax_top.spines['top'].set_visible(False)
+    ax_top.spines['right'].set_visible(False)
+    ax_top.spines['left'].set_visible(False)
+    ax_top.grid(False)
+
+    # Right marginal
+    ax_right.tick_params(axis='y', labelleft=False, left=False)
+    ax_right.tick_params(axis='x', bottom=False, labelbottom=False)
+    ax_right.spines['top'].set_visible(False)
+    ax_right.spines['right'].set_visible(False)
+    ax_right.spines['bottom'].set_visible(False)
+    ax_right.grid(False)
+
+def add_residual_inset(ax, residuals, color):
+    # Moved upward to avoid overlapping with main x-axis labels
+    ax_inset = ax.inset_axes([0.60, 0.20, 0.28, 0.18])
+
+    ax_inset.hist(residuals, bins=28, color=color, alpha=0.70, edgecolor='gray', linewidth=0.6)
+
+    mean_res = np.mean(residuals)
+    ax_inset.axvline(mean_res, color='k', linestyle='--', linewidth=1.0)
+
+    x_left, x_right = ax_inset.get_xlim()
+    y_bottom, y_top = ax_inset.get_ylim()
+    x_shift = 0.04 * (x_right - x_left)
+
+    ax_inset.text(mean_res + x_shift, y_top * 0.82, f'{mean_res:.2f}',
+                  color=color, fontsize=9, fontweight='bold')
+
+    ax_inset.set_xlabel('Residual', fontsize=9, labelpad=1)
+    ax_inset.set_ylabel('')
+    ax_inset.set_yticks([])
+    ax_inset.tick_params(axis='x', labelsize=8, pad=1)
+
+    ax_inset.spines['top'].set_visible(False)
+    ax_inset.spines['right'].set_visible(False)
+    ax_inset.spines['left'].set_visible(False)
+    ax_inset.patch.set_alpha(0.0)
+
+def plot_parity_with_marginals(ax, true, pred, xlabel, ylabel, color, mae_text, rmse_text):
+    true = np.asarray(true).reshape(-1)
+    pred = np.asarray(pred).reshape(-1)
+
+    xmin, xmax = get_limits(true, pred, padding=0.05)
+
+    ax.scatter(true, pred, s=35, c=color, alpha=0.35, edgecolors='none', rasterized=True)
+    ax.plot([xmin, xmax], [xmin, xmax], color='grey', linestyle='--', linewidth=2.0)
+
+    ax.set_xlim(xmin, xmax)
+    ax.set_ylim(xmin, xmax)
+
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel(ylabel)
+
+    beautify_axes(ax)
+
+    r2_val = r2_score_np(true, pred)
+    ax.text(0.05, 0.96, f'$R^2 = {r2_val:.4f}$',
+            transform=ax.transAxes, color='black', fontsize=11, va='top')
+    ax.text(0.05, 0.86, mae_text,
+            transform=ax.transAxes, color='black', fontsize=11, va='top')
+    ax.text(0.05, 0.76, rmse_text,
+            transform=ax.transAxes, color='black', fontsize=11, va='top')
+
+    add_marginal_distributions(ax, true, pred, color, bins=35)
+
+    residuals = pred - true
+    add_residual_inset(ax, residuals, color)
+
+# =========================
+# Figure and colors - Changed to 1x3 layout
+# =========================
+fig, axs = plt.subplots(1, 3, figsize=(12, 4.0), dpi=100)
+
+energy_color = '#1f77b4'   # blue
+force_color = '#2ca02c'    # green
+stress_color = '#ff7f0e'   # orange
+
+# =========================
+# Energy panel
+# energy_data[:, 0] -> NEP
+# energy_data[:, 1] -> DFT
+# =========================
+energy_true = energy_data[:, 1]
+energy_pred = energy_data[:, 0]
+energy_rmse = rmse(energy_pred, energy_true) * 1000.0
+energy_mae = mae(energy_pred, energy_true) * 1000.0
+
+plot_parity_with_marginals(
+    axs[0],
+    energy_true,
+    energy_pred,
+    'DFT energy (eV/atom)',
+    'NEP energy (eV/atom)',
+    energy_color,
+    f'MAE = {energy_mae:.2f} meV/atom', 
+    f'RMSE = {energy_rmse:.2f} meV/atom'
+)
+
+# =========================
+# Force panel
+# force_data[:, 0:3] -> NEP
+# force_data[:, 3:6] -> DFT
+# =========================
+force_true = force_data[:, 3:6]
+force_pred = force_data[:, 0:3]
+
+force_true_flat = force_true.reshape(-1)
+force_pred_flat = force_pred.reshape(-1)
+
+force_rmse = rmse(force_pred_flat, force_true_flat) * 1000.0
+force_mae = mae(force_pred_flat, force_true_flat) * 1000.0
+
+plot_parity_with_marginals(
+    axs[1],
+    force_true_flat,
+    force_pred_flat,
+    r'DFT force (eV/$\mathrm{\AA}$)',
+    r'NEP force (eV/$\mathrm{\AA}$)',
+    force_color,
+    rf'MAE = {force_mae:.2f} meV/$\mathrm{{\AA}}$',
+    rf'RMSE = {force_rmse:.2f} meV/$\mathrm{{\AA}}$'
+)
+
+# =========================
+# Stress panel
+# stress_data[:, 0:6] -> NEP
+# stress_data[:, 6:12] -> DFT
+# =========================
+stress_true = stress_data[:, 6:12]
+stress_pred = stress_data[:, 0:6]
+
+stress_true_flat = stress_true.reshape(-1)
+stress_pred_flat = stress_pred.reshape(-1)
+
+stress_rmse = rmse(stress_pred_flat, stress_true_flat)
+stress_mae = mae(stress_pred_flat, stress_true_flat)
+
+plot_parity_with_marginals(
+    axs[2],
+    stress_true_flat,
+    stress_pred_flat,
+    'DFT stress (GPa)',
+    'NEP stress (GPa)',
+    stress_color,
+    f'MAE = {stress_mae:.4f} GPa',
+    f'RMSE = {stress_rmse:.4f} GPa'
+)
+
+# =========================
+# Layout and output
+# =========================
+abc = [f"({chr(97 + i)})" for i in range(3)]
+
+for i, ax in enumerate(axs.flat):
+    ax.text(-0.15, 1.05, abc[i],
+            transform=ax.transAxes,
+            fontsize=15,
+            ha='left',
+            va='bottom')
+
 plt.tight_layout()
-#fig.subplots_adjust(top=0.968,bottom=0.16, left=0.086,right=0.983,hspace=0.2,wspace=0.25)
+plt.subplots_adjust(wspace=0.2)
 
 if len(sys.argv) > 1 and sys.argv[1] == 'save':
-    plt.savefig('prediction.png', dpi=300)
+    plt.savefig('prediction.png', dpi=300, bbox_inches='tight')
 else:
-    # Check if the current backend is non-interactive
     from matplotlib import get_backend
     if get_backend().lower() in ['agg', 'cairo', 'pdf', 'ps', 'svg']:
-        print("Unable to display the plot due to the non-interactive backend.")
-        print("The plot has been automatically saved as 'prediction.png'.")
-        plt.savefig('prediction.png', dpi=300)
+        print("Non-interactive backend detected. Plot saved as 'prediction.png'.")
+        plt.savefig('prediction.png', dpi=300, bbox_inches='tight')
     else:
         plt.show()
diff --git a/Scripts/plt_scripts/plt_prediction2.py b/Scripts/plt_scripts/plt_prediction2.py
new file mode 100644
index 0000000..aa27ae2
--- /dev/null
+++ b/Scripts/plt_scripts/plt_prediction2.py
@@ -0,0 +1,104 @@
+import sys
+import matplotlib.pyplot as plt
+import numpy as np
+# from cycler import cycler
+
+# custom_colors = ['#8FA2CD', '#F8BC7E', '#A9CA70', '#F09BA0', '#9D9EA3', '#B7B7EB']
+# plt.rcParams['axes.prop_cycle'] = cycler(color=custom_colors)
+
+# Load data
+energy_data = np.loadtxt('energy_train.out')
+force_data = np.loadtxt('force_train.out')
+#virial_data = np.loadtxt('virial_train.out')
+stress_data = np.loadtxt('stress_train.out')
+
+# Filter out rows with invalid stress data
+valid_rows = ~np.any(np.abs(stress_data[:, :12]) >= 1e6, axis=1)
+stress_data = stress_data[valid_rows]
+
+# Function to calculate RMSE
+def calculate_rmse(pred, actual):
+    return np.sqrt(np.mean((pred - actual) ** 2))
+
+# Function to calculate dynamic axis limits
+def calculate_limits(train_data, padding=0.08):
+    data_min = np.min(train_data)
+    data_max = np.max(train_data)
+    data_range = data_max - data_min
+    return data_min - padding * data_range, data_max + padding * data_range
+
+# Create a subplot with 1 row and 3 columns
+fig, axs = plt.subplots(1, 3, figsize=(12, 3.3), dpi=100)
+
+# Plotting the train_data figure
+xmin_energy, xmax_energy = calculate_limits(energy_data[:, 1])
+axs[0].set_xlim(xmin_energy, xmax_energy)
+axs[0].set_ylim(xmin_energy, xmax_energy)
+axs[0].plot(energy_data[:, 1], energy_data[:, 0], '.', markersize=10)
+axs[0].plot([xmin_energy, xmax_energy], [xmin_energy, xmax_energy], linewidth=2, color='grey', linestyle='--')
+axs[0].set_xlabel('DFT energy (eV/atom)', fontsize=11)
+axs[0].set_ylabel('NEP energy (eV/atom)', fontsize=11)
+axs[0].legend(['energy'], loc='upper left')
+axs[0].axis('tight')
+axs[0].tick_params(axis='both', labelsize=11)
+
+# Calculate and display RMSE for energy
+energy_rmse = calculate_rmse(energy_data[:, 0], energy_data[:, 1]) * 1000
+axs[0].text(0.35, 0.08, f'RMSE: {energy_rmse:.2f} meV/atom', transform=axs[0].transAxes, fontsize=11, verticalalignment='center')
+#axs[0].text(-0.1, 1.03, "(a)", transform=axs[0].transAxes, fontsize=13, va='top', ha='right')
+
+# Plotting the force_data figure
+xmin_force, xmax_force = calculate_limits(force_data[:, 3:6].reshape(-1))
+axs[1].set_xlim(xmin_force, xmax_force)
+axs[1].set_ylim(xmin_force, xmax_force)
+axs[1].plot(force_data[:, 3:6], force_data[:, 0:3], '.', markersize=10)
+axs[1].plot([xmin_force, xmax_force], [xmin_force, xmax_force], linewidth=2, color='grey', linestyle='--')
+axs[1].set_xlabel(r'DFT force (eV/$\mathrm{\AA}$)', fontsize=11)
+axs[1].set_ylabel(r'NEP force (eV/$\mathrm{\AA}$)', fontsize=11)
+axs[1].tick_params(axis='both', labelsize=11)
+axs[1].legend(['fx', 'fy', 'fz'], loc='upper left')
+axs[1].axis('tight')
+
+# Calculate and display RMSE for forces
+force_pred_flat = force_data[:, 0:3].flatten()
+force_actual_flat = force_data[:, 3:6].flatten()
+force_rmse = calculate_rmse(force_pred_flat, force_actual_flat) * 1000
+mean_force_rmse = force_rmse
+axs[1].text(0.35, 0.08, rf'RMSE: {mean_force_rmse:.2f} meV/$\mathrm{{\AA}}$', transform=axs[1].transAxes, fontsize=11, verticalalignment='center')
+#axs[1].text(-0.1, 1.03, "(b)", transform=axs[1].transAxes, fontsize=13, va='top', ha='right')
+
+# Plotting the stress figure
+xmin_stress, xmax_stress = calculate_limits(stress_data[:, 6:12].reshape(-1))
+axs[2].set_xlim(xmin_stress, xmax_stress)
+axs[2].set_ylim(xmin_stress, xmax_stress)
+axs[2].plot(stress_data[:, 6:12], stress_data[:, 0:6], '.', markersize=10)
+axs[2].plot([xmin_stress, xmax_stress], [xmin_stress, xmax_stress], linewidth=2, color='grey', linestyle='--')
+axs[2].set_xlabel('DFT stress (GPa)', fontsize=11)
+axs[2].set_ylabel('NEP stress (GPa)', fontsize=11)
+axs[2].tick_params(axis='both', labelsize=11)
+axs[2].legend(['xx', 'yy', 'zz', 'xy', 'yz', 'zx'], loc='upper left')
+axs[2].axis('tight')
+
+# Calculate and display RMSE for stresses
+stress_pred_flat = stress_data[:, 0:6].flatten()
+stress_actual_flat = stress_data[:, 6:12].flatten()
+stress_rmse = calculate_rmse(stress_pred_flat, stress_actual_flat)
+mean_stress_rmse = stress_rmse
+axs[2].text(0.35, 0.08, f'RMSE: {mean_stress_rmse:.4f} GPa', transform=axs[2].transAxes, fontsize=11, verticalalignment='center')
+#axs[2].text(-0.1, 1.03, "(c)", transform=axs[2].transAxes, fontsize=13, va='top', ha='right')
+
+# Adjust layout for better spacing
+plt.tight_layout()
+#fig.subplots_adjust(top=0.968,bottom=0.16, left=0.086,right=0.983,hspace=0.2,wspace=0.25)
+
+if len(sys.argv) > 1 and sys.argv[1] == 'save':
+    plt.savefig('prediction.png', dpi=300)
+else:
+    # Check if the current backend is non-interactive
+    from matplotlib import get_backend
+    if get_backend().lower() in ['agg', 'cairo', 'pdf', 'ps', 'svg']:
+        print("Unable to display the plot due to the non-interactive backend.")
+        print("The plot has been automatically saved as 'prediction.png'.")
+        plt.savefig('prediction.png', dpi=300)
+    else:
+        plt.show()
diff --git a/Scripts/plt_scripts/plt_train_density.py b/Scripts/plt_scripts/plt_train_density.py
index af92b93..0a37f35 100644
--- a/Scripts/plt_scripts/plt_train_density.py
+++ b/Scripts/plt_scripts/plt_train_density.py
@@ -61,7 +61,7 @@ def calculate_limits(train_data, padding=0.08):
 axs[0, 1].set_xlim(xmin_energy, xmax_energy)
 axs[0, 1].set_ylim(xmin_energy, xmax_energy)
 
-axs[0, 1].hist2d(energy_data[:, 0], energy_data[:, 1], 
+axs[0, 1].hist2d(energy_data[:, 1], energy_data[:, 0], 
                  bins=100, cmap='Blues', cmin=1, norm=LogNorm(), 
                  range=[[xmin_energy, xmax_energy], [xmin_energy, xmax_energy]])
 
@@ -77,8 +77,8 @@ def calculate_limits(train_data, padding=0.08):
                transform=axs[0, 1].transAxes, fontsize=10, verticalalignment='center', horizontalalignment='center')
 
 
-f_target = force_data[:, 0:3].flatten()
-f_pred = force_data[:, 3:6].flatten()
+f_pred = force_data[:, 0:3].flatten()
+f_target = force_data[:, 3:6].flatten()
 
 xmin_force, xmax_force = calculate_limits(f_target)
 axs[1, 0].set_xlim(xmin_force, xmax_force)
@@ -105,8 +105,8 @@ def calculate_limits(train_data, padding=0.08):
 if stress_data.shape[0] == 0:
     axs[1, 1].axis('off')
 else:
-    s_target = stress_data[:, 0:6].flatten()
-    s_pred = stress_data[:, 6:12].flatten()
+    s_pred = stress_data[:, 0:6].flatten()
+    s_target = stress_data[:, 6:12].flatten()
 
     xmin_stress, xmax_stress = calculate_limits(s_target)
     axs[1, 1].set_xlim(xmin_stress, xmax_stress)
diff --git a/Scripts/sample_structures/neptrain_select_structs.py b/Scripts/sample_structures/neptrain_select_structs.py
index a4b738a..973506e 100644
--- a/Scripts/sample_structures/neptrain_select_structs.py
+++ b/Scripts/sample_structures/neptrain_select_structs.py
@@ -148,6 +148,10 @@ def FarthestPointSample(new_data, now_data=[], min_distance=None, min_select=1,
     side_kde.set_xticks([])
     side_kde.set_yticks([])
 
+np.savetxt('pca_sample.txt', proj_sample, fmt='%.8f', header='sample_x sample_y', comments='')
+np.savetxt('pca_train.txt', proj_train, fmt='%.8f', header='train_x train_y', comments='')
+np.savetxt('pca_selected.txt', proj_selected, fmt='%.8f', header='selected_x selected_y', comments='')
+
 plt.tight_layout()
 #plt.show()
 plt.savefig('select.png')
diff --git a/docs/Gallery/prediction.png b/docs/Gallery/prediction.png
index 95d8aa2..fd4d2fe 100644
Binary files a/docs/Gallery/prediction.png and b/docs/Gallery/prediction.png differ
diff --git a/docs/publications.json b/docs/publications.json
index cdd95f3..ef50505 100644
--- a/docs/publications.json
+++ b/docs/publications.json
@@ -1,4 +1,11 @@
 [
+    {
+        "title": "Bridging accuracy and efficiency in solid-liquid interface simulations: An accelerated fine-tuning strategy for MgO hydration",
+        "authors": "Yuanrui Wang, Sijin Zhou, Ziliang Wang",
+        "journal": "Comput. Mater. Sci. 2026, 269, 114726",
+        "year": "2026",
+        "link": "https://doi.org/10.1016/j.commatsci.2026.114726"
+    },
     {
         "title": "A Perspective on Training Machine Learning Force Fields for Solid-State Electrolyte Materials",
         "authors": "Zihan Yan, Shengjie Tang, Yizhou Zhu",
diff --git a/docs/tutorials/analyzer.md b/docs/tutorials/analyzer.md
index a7e0e1a..ec4480c 100644
--- a/docs/tutorials/analyzer.md
+++ b/docs/tutorials/analyzer.md
@@ -26,11 +26,21 @@ gpumdkit.sh
 You'll see the following menu:
 
 ```
- ------------>>
- 501) Analyze composition of extxyz
- 502) Find outliers of extxyz
- 000) Return to the main menu
- ------------>>
+ +------------------------------------------------------+
+ |                    ANALYZER TOOLS                    |
+ +------------------------------------------------------+
+ | 501) Analyze composition of extxyz                   |
+ | 502) Find outliers of extxyz                         |
+ | 503) Analyze chemical species of extxyz              |
+ | 504) Check charge balance of extxyz                  |
+ | 505) Analyze energy/force/virial range               |
+ | 506) Filter structures by minimum distance           |
+ | 507) Get minimum interatomic distance                |
+ | 508) Probability density analysis                    |
+ +------------------------------------------------------+
+ | 000) Return to the main menu                         |
+ +------------------------------------------------------+
+ Input the function number:
 ```
 
 for `501`:
@@ -61,6 +71,84 @@ for `502`:
  Enter stress RMSE threshold (GPa): 0.03
 ```
 
+for `503`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in analyzer      |
+ | Script: analyze_chem_species.py                 |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <input.xyz> you want to analyze
+ Example: train.xyz
+ ------------>>
+```
+
+for `504`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in analyzer      |
+ | Script: charge_balance_check.py                 |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <input.xyz> you want to analyze
+ Examp: train.xyz
+ ------------>>
+```
+
+for `505`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in analyzer      |
+ | Script: energy_force_virial_analyzer.py         |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <input.xyz> and property (energy/force/virial)
+ Example: train.xyz force
+ ------------>>
+```
+
+for `506`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in analyzer      |
+ | Script: filter_structures_by_distance.py        |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <input.xyz> and <min_dist>
+ Example: dump.xyz 1.5
+ ------------>>
+```
+
+for `507`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in analyzer      |
+ | Script: get_min_dist.py / get_min_dist_pbc.py   |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <input.xyz> you want to analyze
+ Example: dump.xyz
+ ------------>>
+```
+
+for `508`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in analyzer      |
+ | Script: probability_density_analysis.py         |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <ref_struct> <trajectory_file> <species> <interval>
+ Example: LLZO.vasp dump.xyz Li 0.25
+ ------------>>
+```
+
 Follow the prompts to complete the function.
 
 ## Command-Line Mode
diff --git a/docs/tutorials/calculators.md b/docs/tutorials/calculators.md
index 9aadebd..c47302c 100644
--- a/docs/tutorials/calculators.md
+++ b/docs/tutorials/calculators.md
@@ -26,14 +26,24 @@ gpumdkit.sh
 You'll see the following menu:
 
 ```
- ------------>>
- 401) Calc ionic conductivity
- 402) Calc properties by nep
- 403) Calc descriptors of specific elements
- 404) Calc density of atomistic states (DOAS)
- 405) Calc nudged elastic band (NEB) by nep
- 000) Return to the main menu
- ------------>>
+ +----------------------------------------------------------+
+ |                     CALCULATOR TOOLS                     |
+ +----------------------------------------------------------+
+ | 401) Calc ionic conductivity                             |
+ | 402) Calc properties by nep                              |
+ | 403) Calc descriptors of specific elements               |
+ | 404) Calc density of atomistic states (DOAS)             |
+ | 405) Calc nudged elastic band (NEB) by nep               |
+ | 406) Build neighbor list                                 |
+ | 407) Calc displacement from trajectory                   |
+ | 408) Calc averaged structure                             |
+ | 409) Calc octahedral tilt                                |
+ | 410) Calc polarization for ABO3                          |
+ | 411) Minimize structure by nep                           |
+ | 412) Calc mean square displacement (MSD) from trajectory |
+ +----------------------------------------------------------+
+ | 000) Return to the main menu                             |
+ +----------------------------------------------------------+
  Input the function number:
 ```
 
@@ -101,6 +111,41 @@ for `405`:
  ------------>>
 ```
 
+for `406-410`:
+
+See https://gpumdkit.cn/htmls/polar_material_analysis.html for details
+
+for `411`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in calculators   |
+ | Script: minimize_structure_by_nep.py            |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <structure_file> <nep.txt> [fmax=0.01] [max_steps=1000]
+ Supported structure formats: POSCAR, .xyz
+ Optional arguments:
+  fmax: Force convergence threshold in eV/Ang (default: 0.01)
+  max_steps: Maximum number of optimization steps (default: 1000)
+ Example: POSCAR nep.txt 0.01 1000
+ ------------>>
+```
+
+for `412`:
+
+```
+ >-------------------------------------------------<
+ | This function calls the script in calculators   |
+ | Script: calc_msd.py                             |
+ | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
+ >-------------------------------------------------<
+ Input <extxyz_file> <element_symbol> <dt_fs> [max_corr_steps]
+   Optional argument: max_corr_steps (default: frame number)
+ Example: dump.xyz Li 10
+ ------------>>
+```
+
 Follow the prompts to complete the function.
 
 ## Command-Line Mode
diff --git a/docs/tutorials/format_conversion.md b/docs/tutorials/format_conversion.md
index ae1e471..5795a79 100644
--- a/docs/tutorials/format_conversion.md
+++ b/docs/tutorials/format_conversion.md
@@ -17,15 +17,23 @@ gpumdkit.sh
 You'll see the following menu:
 
 ```
- ------------>>
- 101) Convert VASP to extxyz
- 102) Convert mtp to extxyz
- 103) Convert CP2K to extxyz
- 104) Convert ABACUS to extxyz
- 105) Convert extxyz to POSCAR
- 000) Return to the main menu
- ------------>>
- Input the function number:
+ +-------------------------------------------------------------+
+ |                   FORMAT CONVERSION TOOLS                   |
+ +-------------------------------------------------------------+
+ | 101) VASP to extxyz            106) Add group labels        |
+ | 102) MTP to extxyz             107) Add weight to extxyz    |
+ | 103) CP2K to extxyz            108) Extract frame extxyz    |
+ | 104) ABACUS to extxyz          109) Clean XYZ info          |
+ | 105) extxyz to POSCAR          110) Replicate structure     |
+ +-------------------------------------------------------------+
+ | out2exyz) OUTCAR to extxyz     xdat2exyz) XDATCAR to extxyz |
+ | pos2exyz) POSCAR to extxyz     pos2lmp)   POSCAR to LAMMPS  |
+ | cif2pos)  CIF to POSCAR        lmp2exyz)  LAMMPS to extxyz  |
+ | cif2exyz) CIF to extxyz        traj2exyz) ASE traj to extxyz|
+ +-------------------------------------------------------------+
+ | 000) Return to main menu                                    |
+ +-------------------------------------------------------------+
+ Input the function number or converter keyword:
 ```
 
 for `101`:
diff --git a/docs/tutorials/plot_scripts.md b/docs/tutorials/plot_scripts.md
index 8f5a78b..50d8228 100644
--- a/docs/tutorials/plot_scripts.md
+++ b/docs/tutorials/plot_scripts.md
@@ -15,28 +15,45 @@ gpumdkit.sh -plt <plot_type> [options]
 
 For a complete list:
 ```bash
-gpumdkit.sh -plt -h
-
-+=====================================================================================================+
-|                              GPUMDkit 1.4.2 (dev) (2025-12-17) Plotting Usage                       |
-+=============================================== Plot Types ==========================================+
-| thermo          Plot thermo info                   | train          Plot NEP train results          |
-| prediction      Plot NEP prediction results        | train_test     Plot NEP train and test results |
-| msd             Plot mean square displacement      | msd_conv       Plot the convergence of MSD     |
-| msd_all         Plot MSD of all species            | sdc            Plot self diffusion coefficient |
-| rdf             Plot radial distribution function  | vac            Plot velocity autocorrelation   |
-| restart         Plot parameters in nep.restart     | dimer          Plot dimer plot                 |
-| force_errors    Plot force errors                  | des            Plot descriptors                |
-| charge          Plot charge distribution           | lr             Plot learning rate              |
-| doas            Plot density of atomistic states   | arrhenius_d    Plot Arrhenius diffusivity      |
-| arrhenius_sigma Plot Arrhenius sigma               | net_force      Plot net force distribution     |
-| emd             Plot EMD results                   | nemd           Plot NEMD results               |
-| hnemd           Plot HNEMD results                 |                                                |
-+=====================================================================================================+
-| For detailed usage and examples, use: gpumdkit.sh -plt <plot_type> -h                               |
-+=====================================================================================================+
- See the codes in plt_scripts for more details
- Code path: /d/Westlake/GPUMD/Gpumdkit/Scripts/plt_scripts
+ +-----------------------------------------------------------------------------------------------+
+ |                     GPUMDkit 1.5.5 (dev) (2026-05-10) PLOT & VISUALIZATION TOOLS              |
+ +-----------------------------------------------------------------------------------------------+
+ |  Usage: gpumdkit.sh -plt <type>                        Help: gpumdkit.sh -plt <type> -h       |
+ +-----------------------------------------------------------------------------------------------+
+ |                                    NEP Training & Evaluation                                  |
+ +-----------------------------------------------------------------------------------------------+
+ |  train          - NEP training results           prediction     - NEP prediction results      |
+ |  train_test     - NEP train and test results     parity_density - Parity density plot         |
+ |  train_density  - Training results density plot  restart        - Parameters in nep.restart   |
+ |  charge         - Charge distribution            born_charge    - Born effective charges      |
+ |  dimer          - Dimer energy/force curve       force_errors   - Force errors                |
+ |  des            - Descriptors                    lr             - Learning rate for gnep      |
+ +-----------------------------------------------------------------------------------------------+
+ |                                     Diffusion & Transport                                     |
+ +-----------------------------------------------------------------------------------------------+
+ |  msd            - Mean square displacement       msd_conv       - MSD convergence             |
+ |  msd_all        - MSD for all species            sdc            - Self diffusion coefficient  |
+ |  msd_sdc        - MSD and SDC together           sigma          - Arrhenius ionic conductivity|
+ |  D              - Arrhenius diffusivity          sigma_xyz      - Directional Arrhenius sigma |
+ |  D_xyz          - Directional Arrhenius D                                                     |
+ +-----------------------------------------------------------------------------------------------+
+ |                                    MD & Structural Analysis                                   |
+ +-----------------------------------------------------------------------------------------------+
+ |  thermo         - thermo info in thermo.out      thermo2/3      - Thermo in different styles  |
+ |  rdf            - Radial distribution function   rdf_pmf        - Potential of mean force     |
+ |  vac            - Velocity autocorrelation       cohesive       - Cohesive energy curve       |
+ |  net_force      - Net force distribution         plane-grid     - Displacement plane grid     |
+ |  doas           - Density of atomistic states                                                 |
+ +-----------------------------------------------------------------------------------------------+
+ |                                        Heat Transport                                         |
+ +-----------------------------------------------------------------------------------------------+
+ |  emd            - EMD results                    nemd           - NEMD results                |
+ |  hnemd          - HNEMD results                  viscosity      - Viscosity                   |
+ +-----------------------------------------------------------------------------------------------+
+ |                                          Phonons                                              |
+ +-----------------------------------------------------------------------------------------------+
+ |  pdos           - VAC and PDOS                                                                |
+ +-----------------------------------------------------------------------------------------------+
 ```
 
 ---
diff --git a/docs/tutorials/sample_structures.md b/docs/tutorials/sample_structures.md
index badc110..3778c41 100644
--- a/docs/tutorials/sample_structures.md
+++ b/docs/tutorials/sample_structures.md
@@ -25,14 +25,17 @@ gpumdkit.sh
 You'll see the following menu:
 
 ```
- ------------>>
- 201) Sample structures from extxyz
- 202) Sample structures by pynep
- 203) Sample structures by neptrain
- 204) Perturb structure
- 205) Select max force deviation structs
- 000) Return to the main menu
- ------------>>
+ +------------------------------------------------------+
+ |                 SAMPLE STRUCTURE TOOLS               |
+ +------------------------------------------------------+
+ | 201) Sample structures from extxyz                   |
+ | 202) FPS sampling by PyNEP [deprecated]              |
+ | 203) FPS sampling by NepTrain [preferred]            |
+ | 204) Perturb structure                               |
+ | 205) Select max force deviation structs              |
+ +------------------------------------------------------+
+ | 000) Return to the main menu                         |
+ +------------------------------------------------------+
  Input the function number:
 ```
 
@@ -69,7 +72,7 @@ train.xyz uniform 50
 
 The script `sample_structures.py` in the `Scripts/sample_strcutures` will be called to perform the sampling.
 
-### Option 202: Sample structures by pynep
+### Option 202: FPS sampling by PyNEP
 
 This function calls the `pynep_select_structs.py` in the `Scripts/sample_structures` to sampling the structures by `pynep`.
 
@@ -83,38 +86,17 @@ You will see the following prompt:
 
 ```
  +-------------------------------------------------+
- | PyNEP package is no longer actively maintained  |
- |       Recommend using function 203 instead      |
- +-------------------------------------------------+
+ |  Function 202 (pynep sampling) is deprecated.   |
+ |  PyNEP package is no longer actively maintained.|
+ |  Please use 203) neptrain sampling instead.     |
+ |  If you still insist on pynep, run directly:    |
+ |                 gpumdkit.sh -pynep              |
  +-------------------------------------------------+
- |     To use parallel version, please use:        |
- |            gpumdkit.sh -pynep                   |
- +-------------------------------------------------+
- >-------------------------------------------------<
- | Calling the script in Scripts/sample_structures |
- | Script: pynep_select_structs.py                 |
- | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
- >-------------------------------------------------<
- Input <sample.xyz> <train.xyz> <nep_model>
- Examp: dump.xyz train.xyz nep.txt
- ------------>>
-```
-
-`<samle.xyz>`: extxyz file
-
-`<train.xyz>`: train.xyz
-
-`<nep_model>`: nep.txt
-
-Enter the following parameters:
-
-```
-dump.xyz train.xyz nep.txt
 ```
 
 **Note:** pynep package is no longer actively maintained, this function will be removed soon.
 
-### Option 203: Sample structures by neptrain
+### Option 203: FPS sampling by NepTrain
 
 This function calls the `neptrain_select_structs.py` script in the `Scripts/sample_structures` to sampling the structures by `neptrain`
 
@@ -137,7 +119,17 @@ You will see the following prompt:
  ------------>>
 ```
 
-It is the same with option `202`.
+`<samle.xyz>`: extxyz file
+
+`<train.xyz>`: train.xyz
+
+`<nep_model>`: nep.txt
+
+Enter the following parameters:
+
+```
+dump.xyz train.xyz nep.txt
+```
 
 ### Option 204: Perturb structure
 
@@ -181,7 +173,7 @@ POSCAR 20 0.03 0.2 uniform
 
 The script `perturb_structure.py` in the `Scripts/sample_strcutures` will be called to perform the perturbation.
 
-### Option 205: Select max force deviation structs from active.xyz
+### Option 205: Select max force deviation structs
 
 This option allows you to select max force deviation structures from `active.xyz` generated by the command `active` in `gpumd`.
 
diff --git a/docs/tutorials/tutorials.md b/docs/tutorials/tutorials.md
index 089046b..9a07ad3 100644
--- a/docs/tutorials/tutorials.md
+++ b/docs/tutorials/tutorials.md
@@ -46,20 +46,22 @@ gpumdkit.sh
 
 You'll see:
 ```
-         ____ ____  _   _ __  __ ____  _    _ _
-        / ___|  _ \| | | |  \/  |  _ \| | _(_) |_
-       | |  _| |_) | | | | |\/| | | | | |/ / | __|
-       | |_| |  __/| |_| | |  | | |_| |   <| | |_
-        \____|_|    \___/|_|  |_|____/|_|\_\_|\__|
+           ____ ____  _   _ __  __ ____  _    _ _
+          / ___|  _ \| | | |  \/  |  _ \| | _(_) |_
+         | |  _| |_) | | | | |\/| | | | | |/ / | __|
+         | |_| |  __/| |_| | |  | | |_| |   <| | |_
+          \____|_|    \___/|_|  |_|____/|_|\_\_|\__|
 
-        GPUMDkit Version 1.4.2 (dev) (2025-12-17)
-  Core Developer: Zihan YAN (yanzihan@westlake.edu.cn)
+          GPUMDkit Version 1.5.5 (dev) (2026-05-10)
+    Core Developer: Zihan YAN (yanzihan@westlake.edu.cn)
+ Main Contributors: Denan LI, Xin WU, Zhoulin LIU & Chen HUA
 
- ----------------------- GPUMD -----------------------
+ ---------------------- GPUMD ------------------------
  1) Format Conversion          2) Sample Structures
  3) Workflow                   4) Calculators
- 5) Analyzer                   6) Developing ...
- 0) Quit!
+ 5) Analyzer                   6) Visualization
+ 7) Utilities                  8) Developing...
+ 0) Exit
  ------------>>
  Input the function number:
 ```
@@ -80,28 +82,38 @@ gpumdkit.sh -h
 You will see:
 
 ```
-+==================================================================================================+
-|                              GPUMDkit 1.4.2 (dev) (2025-12-17) Usage                             |
-+======================================== Conversions =============================================+
-| -out2xyz       Convert OUTCAR to extxyz       | -pos2exyz     Convert POSCAR to extxyz           |
-| -cif2pos       Convert cif to POSCAR          | -pos2lmp      Convert POSCAR to LAMMPS           |
-| -cif2exyz      Convert cif to extxyz          | -lmp2exyz     Convert LAMMPS-dump to extxyz      |
-| -addgroup      Add group label                | -addweight    Add weight to the struct in extxyz |
-| Developing...                                 | Developing...                                    |
-+========================================= Analysis ===============================================+
-| -range         Print range of energy etc.     | -max_rmse     Get max RMSE from extxyz           |
-| -min_dist      Get min_dist between atoms     | -min_dist_pbc Get min_dist considering PBC       |
-| -filter_box    Filter struct by box limits    | -filter_value Filter struct by value (efs)       |
-| -filter_dist   Filter struct by min_dist      | -analyze_comp Analyze composition of extxyz      |
-| -pynep         Sample struct by pynep         | Developing...                                    |
-+====================================== Misc Utilities ============================================+
-| -plt           Plot scripts                   | -get_frame     Extract the specified frame       |
-| -calc          Calculators                    | -clean_xyz     Clean extra info in XYZ file      |
-| -clean         Clear files for work_dir       | -time          Time consuming Analyzer           |
-| -update        Update GPUMDkit                | Developing...                                    |
-+==================================================================================================+
-| For detailed usage and examples, use: gpumdkit.sh -<option> -h                                   |
-+==================================================================================================+
++-------------------------------------------------------------------------------------------------------+
+|                          GPUMDkit 1.5.5 (dev) (2026-05-10) Command Help                               |
++-------------------------------------------------------------------------------------------------------+
+|                                          MAIN FUNCTIONS                                               |
++-------------------------------------------------------------------------------------------------------+
+| -h            Show this help table            | -plt <type>        Plot and visualization tools       |
+| -calc <type>  Calculator tools                | -time <gpumd|nep>  Time-consuming analyzer            |
+| -update       Update GPUMDkit                 | -clean             Clean extra files in current dir   |
++-------------------------------------------------------------------------------------------------------+
+|                                         FORMAT CONVERSION                                             |
++-------------------------------------------------------------------------------------------------------+
+| -out2xyz      OUTCAR -> extxyz (shell)        | -out2exyz          OUTCAR -> extxyz (python)          |
+| -cp2k2xyz     CP2K log -> xyz                 | -xdat2exyz         XDATCAR -> extxyz                  |
+| -cif2pos      cif -> POSCAR                   | -cif2exyz          cif -> extxyz                      |
+| -pos2exyz     POSCAR -> extxyz                | -exyz2pos          extxyz -> POSCAR                   |
+| -pos2lmp      POSCAR -> LAMMPS data           | -lmp2exyz          LAMMPS dump -> extxyz              |
+| -traj2exyz    ASE traj -> extxyz              | -replicate         Replicate structure                |
+| -addgroup     Add group labels                | -addweight         Add structure weight in extxyz     |
+| -clean_xyz    Clean extra info in extxyz      | -get_frame         Extract specific frame             |
+| -frame_range  Extract frames by range         |                                                       |
++-------------------------------------------------------------------------------------------------------+
+|                                            ANALYSIS                                                   |
++-------------------------------------------------------------------------------------------------------+
+| -range        Energy/force/virial statistics  | -analyze_comp      Analyze composition                |
+| -chem_species Analyze chemical species        | -cbc               Charge balance check               |
+| -min_dist     Min distance (no PBC)           | -min_dist_pbc      Min distance with PBC              |
+| -filter_dist  Filter by min_dist (no PBC)     | -filter_dist_pbc   Filter by min_dist (PBC)           |
+| -pda          Probability density analysis    | -hbond             Hydrogen-bond analysis             |
+| -pynep        FPS sampling by PyNEP           |                                                       |
++-------------------------------------------------------------------------------------------------------+
+| Detailed usage: gpumdkit.sh -<option> -h    Plot details: gpumdkit.sh -plt <type> -h                  |
++-------------------------------------------------------------------------------------------------------+
 ```
 
 **Best for:**
@@ -119,15 +131,23 @@ Convert structure files between different formats:
 #### Interactive mode (option 1)
 
 ```shell
- ------------>>
- 101) Convert VASP to extxyz
- 102) Convert mtp to extxyz
- 103) Convert CP2K to extxyz
- 104) Convert ABACUS to extxyz
- 105) Convert extxyz to POSCAR
- 000) Return to the main menu
- ------------>>
- Input the function number:
+ +-------------------------------------------------------------+
+ |                   FORMAT CONVERSION TOOLS                   |
+ +-------------------------------------------------------------+
+ | 101) VASP to extxyz            106) Add group labels        |
+ | 102) MTP to extxyz             107) Add weight to extxyz    |
+ | 103) CP2K to extxyz            108) Extract frame extxyz    |
+ | 104) ABACUS to extxyz          109) Clean XYZ info          |
+ | 105) extxyz to POSCAR          110) Replicate structure     |
+ +-------------------------------------------------------------+
+ | out2exyz) OUTCAR to extxyz     xdat2exyz) XDATCAR to extxyz |
+ | pos2exyz) POSCAR to extxyz     pos2lmp)   POSCAR to LAMMPS  |
+ | cif2pos)  CIF to POSCAR        lmp2exyz)  LAMMPS to extxyz  |
+ | cif2exyz) CIF to extxyz        traj2exyz) ASE traj to extxyz|
+ +-------------------------------------------------------------+
+ | 000) Return to main menu                                    |
+ +-------------------------------------------------------------+
+ Input the function number or converter keyword:
 ```
 
 #### Command-line mode
@@ -161,14 +181,17 @@ Select diverse structures for training:
 #### Interactive mode (option 2)
 
 ```shell
- ------------>>
- 201) Sample structures from extxyz
- 202) Sample structures by pynep
- 203) Sample structures by neptrain
- 204) Perturb structure
- 205) Select max force deviation structs
- 000) Return to the main menu
- ------------>>
+ +------------------------------------------------------+
+ |                 SAMPLE STRUCTURE TOOLS               |
+ +------------------------------------------------------+
+ | 201) Sample structures from extxyz                   |
+ | 202) FPS sampling by PyNEP [deprecated]              |
+ | 203) FPS sampling by NepTrain [preferred]            |
+ | 204) Perturb structure                               |
+ | 205) Select max force deviation structs              |
+ +------------------------------------------------------+
+ | 000) Return to the main menu                         |
+ +------------------------------------------------------+
  Input the function number:
 ```
 
@@ -181,12 +204,15 @@ Batch processing for high-throughput calculations:
 #### Interactive mode (option 3)
 
 ```
- ------------>>
- 301) SCF batch pretreatment
- 302) MD sample batch pretreatment (gpumd)
- 303) MD sample batch pretreatment (lmp)
- 000) Return to the main menu
- ------------>>
+ +---------------------------------------------------------+
+ |                      WORKFLOW TOOLS                     |
+ +---------------------------------------------------------+
+ | 301) SCF batch pretreatment                             |
+ | 302) MD sample batch pretreatment (gpumd)               |
+ | 303) MD sample batch pretreatment (lmp)                 |
+ +---------------------------------------------------------+
+ | 000) Return to the main menu                            |
+ +---------------------------------------------------------+
  Input the function number:
 ```
 
@@ -199,14 +225,24 @@ Compute material properties:
 #### Interactive mode (option 4)
 
 ```shell
- ------------>>
- 401) Calc ionic conductivity
- 402) Calc properties by nep
- 403) Calc descriptors of specific elements
- 404) Calc density of atomistic states (DOAS)
- 405) Calc nudged elastic band (NEB) by nep
- 000) Return to the main menu
- ------------>>
+ +----------------------------------------------------------+
+ |                     CALCULATOR TOOLS                     |
+ +----------------------------------------------------------+
+ | 401) Calc ionic conductivity                             |
+ | 402) Calc properties by nep                              |
+ | 403) Calc descriptors of specific elements               |
+ | 404) Calc density of atomistic states (DOAS)             |
+ | 405) Calc nudged elastic band (NEB) by nep               |
+ | 406) Build neighbor list                                 |
+ | 407) Calc displacement from trajectory                   |
+ | 408) Calc averaged structure                             |
+ | 409) Calc octahedral tilt                                |
+ | 410) Calc polarization for ABO3                          |
+ | 411) Minimize structure by nep                           |
+ | 412) Calc mean square displacement (MSD) from trajectory |
+ +----------------------------------------------------------+
+ | 000) Return to the main menu                             |
+ +----------------------------------------------------------+
  Input the function number:
 ```
 
@@ -232,11 +268,20 @@ Validate and filter structure files:
 #### Interactive mode (option 5)
 
 ```
- ------------>>
- 501) Analyze composition of extxyz
- 502) Find outliers of extxyz
- 000) Return to the main menu
- ------------>>
+ +------------------------------------------------------+
+ |                    ANALYZER TOOLS                    |
+ +------------------------------------------------------+
+ | 501) Analyze composition of extxyz                   |
+ | 502) Find outliers of extxyz                         |
+ | 503) Analyze chemical species of extxyz              |
+ | 504) Check charge balance of extxyz                  |
+ | 505) Analyze energy/force/virial range               |
+ | 506) Filter structures by minimum distance           |
+ | 507) Get minimum interatomic distance                |
+ | 508) Probability density analysis                    |
+ +------------------------------------------------------+
+ | 000) Return to the main menu                         |
+ +------------------------------------------------------+
  Input the function number:
 ```
 
diff --git a/docs/tutorials/workflow.md b/docs/tutorials/workflow.md
index 6d13c17..314ab4b 100644
--- a/docs/tutorials/workflow.md
+++ b/docs/tutorials/workflow.md
@@ -22,12 +22,16 @@ gpumdkit.sh
 You'll see the following menu:
 
 ```
- ------------>>
- 301) SCF batch pretreatment
- 302) MD sample batch pretreatment (gpumd)
- 303) MD sample batch pretreatment (lmp)
- 000) Return to the main menu
- ------------>>
+ +---------------------------------------------------------+
+ |                      WORKFLOW TOOLS                     |
+ +---------------------------------------------------------+
+ | 301) SCF batch pretreatment                             |
+ | 302) MD sample batch pretreatment (gpumd)               |
+ | 303) MD sample batch pretreatment (lmp)                 |
+ +---------------------------------------------------------+
+ | 000) Return to the main menu                            |
+ +---------------------------------------------------------+
+ Input the function number:
 ```
 
 ---
diff --git a/gpumdkit.sh b/gpumdkit.sh
index a5ffcc8..d350516 100755
--- a/gpumdkit.sh
+++ b/gpumdkit.sh
@@ -1,5 +1,10 @@
 #!/bin/bash
 
+# GPUMDkit - A User-Friendly Toolkit for GPUMD and NEP
+# Project URL: https://github.com/zhyan0603/GPUMDkit
+
+# Description: Main entry point for GPUMDkit interactive and command-line interface
+
 # You need to set the path of GPUMD and GPUMDkit in your ~/.bashrc, for example
 # export GPUMDkit_path=/home/yanzihan/software/GPUMDkit
 
@@ -10,7 +15,7 @@ if [ -z "$GPUMDkit_path" ]; then
     exit 1
 fi
 
-VERSION="1.5.4 (dev) (2026-04-22)"
+VERSION="1.5.5 (dev) (2026-05-10)"
 
 plt_path="${GPUMDkit_path}/Scripts/plt_scripts"
 analyzer_path="${GPUMDkit_path}/Scripts/analyzer"
@@ -23,11 +28,12 @@ sample_path="${GPUMDkit_path}/Scripts/sample_structures"
 #--------------------- main script ----------------------
 # Show the menu
 function menu(){
-echo " ----------------------- GPUMD -----------------------"
+echo " ---------------------- GPUMD ------------------------"
 echo " 1) Format Conversion          2) Sample Structures   "
 echo " 3) Workflow                   4) Calculators         "
-echo " 5) Analyzer                   6) Developing ...      "
-echo " 0) Quit!"
+echo " 5) Analyzer                   6) Visualization       "
+echo " 7) Utilities                  8) Developing...       "
+echo " 0) Exit                                              "
 }
 
 # Function main
@@ -35,12 +41,13 @@ function main(){
     echo " ------------>>"
     echo ' Input the function number:'
     array_choice=(
-        "0" "1" "101" "102" "103" "104" "105"
+        "0" "1" "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
         "2" "201" "202" "203" "204" "205" 
         "3" "301" "302" "303" 
-        "4" "401" "402" "403" "404" "405" "406" "407" "408" "409" "410"
-        "5" "501" "502"
+        "4" "401" "402" "403" "404" "405" "406" "407" "408" "409" "410" "411" "412"
+        "5" "501" "502" "503" "504" "505" "506" "507" "508"
         "6"
+        "7" "701"
     ) 
     read -p " " choice
     while ! echo "${array_choice[@]}" | grep -wq "$choice" 
@@ -64,6 +71,11 @@ function main(){
                 "103") f103_cp2k2xyz ;;
                 "104") f104_abacus2xyz ;;
                 "105") f105_extxyz2poscar ;;
+                "106") f106_add_group_labels ;;
+                "107") f107_add_weight ;;
+                "108") f108_get_frame ;;
+                "109") f109_clean_xyz ;;
+                "110") f110_replicate_structure ;;
             esac ;;
         "2")
             source ${GPUMDkit_path}/src/f2_sample_structures.sh
@@ -102,6 +114,8 @@ function main(){
                 "408") f408_calc_averaged_structure ;;
                 "409") f409_calc_oct_tilt ;;
                 "410") f410_calc_polarization_abo3 ;;
+                "411") f411_minimize_structure_by_nep ;;
+                "412") f412_calc_msd_from_trajectory ;;
             esac ;;           
         "5")
             source ${GPUMDkit_path}/src/f5_analyzers.sh
@@ -109,9 +123,27 @@ function main(){
                 "5") f5_analyzers ;;
                 "501") f501_analyze_composition ;;
                 "502") f502_find_outliers ;;
+                "503") f503_analyze_chem_species ;;
+                "504") f504_charge_balance_check ;;
+                "505") f505_energy_force_virial_analyzer ;;
+                "506") f506_filter_structures_by_distance ;;
+                "507") f507_get_min_dist ;;
+                "508") f508_probability_density_analysis ;;
             esac ;;  
         "6")
-            echo " Developing ..." ;;
+            source ${GPUMDkit_path}/src/f6_plots.sh
+            terminal_width=$(tput cols 2>/dev/null || echo 80)
+            if [ "$terminal_width" -ge 98 ]; then
+                f6_plots_two_column
+            else
+                f6_plots_one_column
+            fi ;;
+        "7")
+            source ${GPUMDkit_path}/src/f7_utilities.sh
+            case $choice in
+                "7") f7_utilities ;;
+                "701") f701_time_consuming_analyzer ;;
+            esac ;;
         *)
             echo " Incorrect Options"
             ;;
@@ -125,55 +157,38 @@ function main(){
 # It will show the usage of each function
 
 function help_info_table(){
-    echo "+==================================================================================================+"
-    echo "|                              GPUMDkit ${VERSION} Usage                             |"
-    echo "+======================================== Conversions =============================================+"
-    echo "| -out2xyz       Convert OUTCAR to extxyz       | -pos2exyz     Convert POSCAR to extxyz           |"
-    echo "| -cif2pos       Convert cif to POSCAR          | -pos2lmp      Convert POSCAR to LAMMPS           |"
-    echo "| -cif2exyz      Convert cif to extxyz          | -lmp2exyz     Convert LAMMPS-dump to extxyz      |"
-    echo "| -addgroup      Add group label                | -addweight    Add weight to the struct in extxyz |"
-    echo "| -cp2k2xyz      Convert CP2K file to extxyz    | -traj2exyz    Convert ASE traj to extxyz         |"
-    echo "| -xdat2exyz     Convert XDATCAR to extxyz      | Developing...                                    |"
-    echo "+========================================= Analysis ===============================================+"
-    echo "| -range         Print range of energy etc.     | -max_rmse     Get max RMSE from extxyz           |"
-    echo "| -min_dist      Get min_dist between atoms     | -min_dist_pbc Get min_dist considering PBC       |"
-    echo "| -filter_box    Filter struct by box limits    | -filter_value Filter struct by value (efs)       |"
-    echo "| -filter_dist   Filter struct by min_dist      | -analyze_comp Analyze composition of extxyz      |"
-    echo "| -pynep         Sample struct by pynep         | Developing...                                    |"
-    echo "+====================================== Misc Utilities ============================================+"
-    echo "| -plt           Plot scripts                   | -get_frame     Extract the specified frame       |"
-    echo "| -calc          Calculators                    | -frame_range   Extract frames by fraction range  |"
-    echo "| -clean         Clear files for work_dir       | -clean_xyz     Clean extra info in XYZ file      |"
-    echo "| -time          Time consuming Analyzer        | -update        Update GPUMDkit                   |"
-    echo "+==================================================================================================+"
-    echo "| For detailed usage and examples, use: gpumdkit.sh -<option> -h                                   |"
-    echo "+==================================================================================================+"
-}
-
-function plot_info_table(){
-    echo "+=====================================================================================================+"
-    echo "|                              GPUMDkit ${VERSION} Plotting Usage                       |"
-    echo "+=============================================== Plot Types ==========================================+"
-    echo "| thermo          Plot thermo info                   | train          Plot NEP train results          |"
-    echo "| prediction      Plot NEP prediction results        | train_test     Plot NEP train and test results |"
-    echo "| msd             Plot mean square displacement      | msd_conv       Plot the convergence of MSD     |"
-    echo "| msd_all         Plot MSD of all species            | sdc            Plot self diffusion coefficient |"
-    echo "| msd_sdc         Plot MSD and SDC together          | cohesive       Plot cohsive energy             |"
-    echo "| rdf             Plot radial distribution function  | vac            Plot velocity autocorrelation   |"
-    echo "| restart         Plot parameters in nep.restart     | dimer          Plot dimer plot                 |"
-    echo "| force_errors    Plot force errors                  | des            Plot descriptors                |"
-    echo "| charge          Plot charge distribution           | lr             Plot learning rate              |"
-    echo "| doas            Plot density of atomistic states   | net_force      Plot net force distribution     |"
-    echo "| sigma           Plot Arrhenius sigma               | D              Plot Arrhenius diffusivity      |"
-    echo "| sigma_xyz       Plot directional Arrhenius sigma   | D_xyz          Plot directional Arrhenius D    |"
-    echo "| emd             Plot EMD results                   | nemd           Plot NEMD results               |"
-    echo "| hnemd           Plot HNEMD results                 | pdos           Plot VAC and PDOS               |"
-    echo "| plane-grid      Plot displacement plane grid       | parity_density Plot parity plot density        |"
-    echo "| rdf_pmf         Plot potential of mean force (PMF) | viscosity      Plot visconsity                 |"
-    echo "| born_charge     Plot Born effective charges        |                                                |"
-    echo "+=====================================================================================================+"
-    echo "| For detailed usage and examples, use: gpumdkit.sh -plt <plot_type> -h                               |"
-    echo "+=====================================================================================================+"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "|                          GPUMDkit ${VERSION} Command Help                               |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "|                                          MAIN FUNCTIONS                                               |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "| -h            Show this help table            | -plt <type>        Plot and visualization tools       |"
+    echo "| -calc <type>  Calculator tools                | -time <gpumd|nep>  Time-consuming analyzer            |"
+    echo "| -update       Update GPUMDkit                 | -clean             Clean extra files in current dir   |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "|                                         FORMAT CONVERSION                                             |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "| -out2xyz      OUTCAR -> extxyz (shell)        | -out2exyz          OUTCAR -> extxyz (python)          |"
+    echo "| -cp2k2xyz     CP2K log -> xyz                 | -xdat2exyz         XDATCAR -> extxyz                  |"
+    echo "| -cif2pos      cif -> POSCAR                   | -cif2exyz          cif -> extxyz                      |"
+    echo "| -pos2exyz     POSCAR -> extxyz                | -exyz2pos          extxyz -> POSCAR                   |"
+    echo "| -pos2lmp      POSCAR -> LAMMPS data           | -lmp2exyz          LAMMPS dump -> extxyz              |"
+    echo "| -traj2exyz    ASE traj -> extxyz              | -replicate         Replicate structure                |"
+    echo "| -addgroup     Add group labels                | -addweight         Add structure weight in extxyz     |"
+    echo "| -clean_xyz    Clean extra info in extxyz      | -get_frame         Extract specific frame             |"
+    echo "| -frame_range  Extract frames by range         |                                                       |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "|                                            ANALYSIS                                                   |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "| -range        Energy/force/virial statistics  | -analyze_comp      Analyze composition                |"
+    echo "| -chem_species Analyze chemical species        | -cbc               Charge balance check               |"
+    echo "| -min_dist     Min distance (no PBC)           | -min_dist_pbc      Min distance with PBC              |"
+    echo "| -filter_dist  Filter by min_dist (no PBC)     | -filter_dist_pbc   Filter by min_dist (PBC)           |"
+    echo "| -pda          Probability density analysis    | -hbond             Hydrogen-bond analysis             |"
+    echo "| -pynep        FPS sampling by PyNEP           |                                                       |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
+    echo "| Detailed usage: gpumdkit.sh -<option> -h    Plot details: gpumdkit.sh -plt <type> -h                  |"
+    echo "+-------------------------------------------------------------------------------------------------------+"
 }
 
 # function citation(){
@@ -270,10 +285,11 @@ if [ ! -z "$1" ]; then
                         echo " | set and charge_train.out.                                |"
                         echo " +----------------------------------------------------------+"
                         python ${plt_path}/plt_charge.py $3 ;;
-                    *) plot_info_table; exit 1 ;;
+                    *) source ${GPUMDkit_path}/src/f6_plots.sh; f6_plots_two_column; exit 1 ;;
                 esac
             else
-                plot_info_table
+                source ${GPUMDkit_path}/src/f6_plots.sh
+                f6_plots_two_column
                 echo " See the codes in plt_scripts for more details"
                 echo " Code path: ${GPUMDkit_path}/Scripts/plt_scripts"
             fi ;;
@@ -329,6 +345,14 @@ if [ ! -z "$1" ]; then
                             echo " Code path: ${calc_path}/calc_doas.py"
                             exit 1
                         fi ;;
+                    minimize)
+                        echo " Calling script by Zihan YAN. "
+                        echo " Code path: ${calc_path}/calc_minimize.py"
+                        python ${calc_path}/calc_minimize.py ${@:3} ;;                     
+                    msd)
+                        echo " Calling script by Zihan YAN. "
+                        echo " Code path: ${calc_path}/calc_msd.py"
+                        python ${calc_path}/calc_msd.py ${@:3} ;;                        
                     nlist)
                         echo " Calling script by Denan LI. "
                         echo " Code path: ${calc_path}/calc_neighbor_list.py"
@@ -459,7 +483,7 @@ if [ ! -z "$1" ]; then
                 echo " Code path: ${format_conv_path}/xdatcar2exyz.py"
                 python ${format_conv_path}/xdatcar2exyz.py $2 $3
             else
-                echo " Usage: -xdat2exyz XDATCAR model.xyz"
+                echo " Usage: -xdat2exyz XDATCAR dump.xyz"
                 echo " See the source code of xdatcar2exyz.py for more details"
                 echo " Code path: ${format_conv_path}/xdatcar2exyz.py"
             fi ;;
diff --git a/src/f1_format_conversions.sh b/src/f1_format_conversions.sh
index e24abae..a407ac4 100644
--- a/src/f1_format_conversions.sh
+++ b/src/f1_format_conversions.sh
@@ -160,24 +160,257 @@ echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/exyz2pos.py"
 echo " ---------------------------------------------------"
 }
 
-#--------------------- function 1 ----------------------
-function f1_format_conversion(){
+# Add group labels
+function f106_add_group_labels(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: add_groups.py                           |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <POSCAR> <element1> <element2> ..."
+echo " Example: POSCAR Li Y Cl"
+echo " ------------>>"
+read -r -a addgroup_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/add_groups.py "${addgroup_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/add_groups.py"
+echo " ---------------------------------------------------"
+}
+
+# Add weight to extxyz
+function f107_add_weight(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: add_weight.py                           |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> <output.xyz> <weight>"
+echo " Example: train.xyz train_weighted.xyz 5"
+echo " ------------>>"
+read -r -a addweight_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/add_weight.py "${addweight_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/add_weight.py"
+echo " ---------------------------------------------------"
+}
+
+# Get a single frame from extxyz
+function f108_get_frame(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: get_frame.py                            |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> <frame_index>"
+echo " Example: train.xyz 1"
+echo " ------------>>"
+read -r -a getframe_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/get_frame.py "${getframe_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/get_frame.py"
+echo " ---------------------------------------------------"
+}
+
+# Clean extra info in XYZ file
+function f109_clean_xyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: clean_xyz.py                            |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> <output.xyz>"
+echo " Example: train.xyz train_clean.xyz"
+echo " ------------>>"
+read -r -a cleanxyz_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/clean_xyz.py "${cleanxyz_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/clean_xyz.py"
+echo " ---------------------------------------------------"
+}
+
+# Replicate structure
+function f110_replicate_structure(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: replicate.py                            |"
+echo " | Developer: Boyi SITU (situboyi@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <inputfile> <outputfile> <a b c / target_num>"
+echo " Example 1: POSCAR POSCAR_222.vasp 2 2 2"
+echo " Example 2: POSCAR POSCAR_256.vasp 256"
+echo " ------------>>"
+read -r -a replicate_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/replicate.py "${replicate_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/replicate.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert OUTCAR to extxyz (Python version)
+function out2exyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: out2exyz.py                             |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input the directory containing OUTCARs"
+echo " Example: ./ "
+echo " ------------>>"
+read -r -p " " dir_outcars
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/out2exyz.py "${dir_outcars}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/out2exyz.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert POSCAR to extxyz
+function pos2extxyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: pos2exyz.py                             |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <POSCAR> <output.xyz>"
+echo " Example: POSCAR model.xyz"
+echo " ------------>>"
+read -r -a pos2extxyz_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/pos2exyz.py "${pos2extxyz_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/pos2exyz.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert CIF to POSCAR
+function cif2poscar(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: cif2pos.py                              |"
+echo " | Developer: Boyi SITU (situboyi@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.cif> <output.vasp>"
+echo " Example: input.cif POSCAR.vasp"
+echo " ------------>>"
+read -r -a cif2pos_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/cif2pos.py "${cif2pos_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/cif2pos.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert CIF to extxyz
+function cif2extxyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: cif2exyz.py                             |"
+echo " | Developer: Boyi SITU (situboyi@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.cif> <output.xyz>"
+echo " Example: input.cif model.xyz"
+echo " ------------>>"
+read -r -a cif2extxyz_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/cif2exyz.py "${cif2extxyz_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/cif2exyz.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert XDATCAR to extxyz
+function xdatcar2extxyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: xdatcar2exyz.py                         |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <XDATCAR> <output.xyz>"
+echo " Example: XDATCAR dump.xyz"
+echo " ------------>>"
+read -r -a xdatcar_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/xdatcar2exyz.py "${xdatcar_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/xdatcar2exyz.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert POSCAR to LAMMPS data
+function poscar2lammps(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: pos2lmp.py                              |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <POSCAR> <lammps.data>"
+echo " Example: POSCAR lammps.data"
 echo " ------------>>"
-echo " 101) Convert VASP to extxyz"
-echo " 102) Convert mtp to extxyz"
-echo " 103) Convert CP2K to extxyz"
-echo " 104) Convert ABACUS to extxyz"
-echo " 105) Convert extxyz to POSCAR"
-echo " 000) Return to the main menu"
+read -r -a pos2lmp_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/pos2lmp.py "${pos2lmp_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/pos2lmp.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert LAMMPS dump to extxyz
+function lammps2extxyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: lmp2exyz.py                             |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <dump_file> <element1> <element2> ..."
+echo " Example: dump.lammpstrj Li O"
 echo " ------------>>"
-echo " Input the function number:"
+read -r -a lmp2extxyz_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/lmp2exyz.py "${lmp2extxyz_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/lmp2exyz.py"
+echo " ---------------------------------------------------"
+}
+
+# Convert ASE traj to extxyz
+function traj2extxyz(){
+echo " >-------------------------------------------------<"
+echo " | Calling the script in Scripts/format_conversion |"
+echo " | Script: traj2exyz.py                            |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.traj> <output.xyz>"
+echo " Example: input.traj output.xyz"
+echo " ------------>>"
+read -r -a traj2extxyz_args
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/format_conversion/traj2exyz.py "${traj2extxyz_args[@]}"
+echo " Code path: ${GPUMDkit_path}/Scripts/format_conversion/traj2exyz.py"
+echo " ---------------------------------------------------"
+}
+
+#--------------------- function 1 ----------------------
+function f1_format_conversion(){
+    echo " +-------------------------------------------------------------+"
+    echo " |                   FORMAT CONVERSION TOOLS                   |"
+    echo " +-------------------------------------------------------------+"
+    echo " | 101) VASP to extxyz            106) Add group labels        |"
+    echo " | 102) MTP to extxyz             107) Add weight to extxyz    |"
+    echo " | 103) CP2K to extxyz            108) Extract frame extxyz    |"
+    echo " | 104) ABACUS to extxyz          109) Clean XYZ info          |"
+    echo " | 105) extxyz to POSCAR          110) Replicate structure     |"
+    echo " +-------------------------------------------------------------+"
+    echo " | out2exyz) OUTCAR to extxyz     xdat2exyz) XDATCAR to extxyz |"
+    echo " | pos2exyz) POSCAR to extxyz     pos2lmp)   POSCAR to LAMMPS  |"
+    echo " | cif2pos)  CIF to POSCAR        lmp2exyz)  LAMMPS to extxyz  |"
+    echo " | cif2exyz) CIF to extxyz        traj2exyz) ASE traj to extxyz|"
+    echo " +-------------------------------------------------------------+"
+    echo " | 000) Return to main menu                                    |"
+    echo " +-------------------------------------------------------------+"
+    echo " Input the function number or converter keyword:"
 
-arry_num_choice=("000" "101" "102" "103" "104" "105" "106") 
+valid_menu_choices=(
+    "000" "101" "102" "103" "104" "105" "106" "107" "108" "109" "110" \
+    "out2exyz" "pos2exyz" "cif2pos" "cif2exyz" "xdat2exyz" "pos2lmp" "lmp2exyz" "traj2exyz"
+) 
 read -p " " num_choice
-while ! echo "${arry_num_choice[@]}" | grep -wq "$num_choice" 
+while ! echo "${valid_menu_choices[@]}" | grep -wq "$num_choice" 
 do
   echo " ------------>>"
-  echo " Please reinput function number..."
+  echo " Please reinput function number or converter keyword..."
   read -p " " num_choice
 done
 
@@ -187,6 +420,19 @@ case $num_choice in
     "103") f103_cp2k2xyz ;;
     "104") f104_abacus2xyz ;;
     "105") f105_extxyz2poscar ;;
+    "106") f106_add_group_labels ;;
+    "107") f107_add_weight ;;
+    "108") f108_get_frame ;;
+    "109") f109_clean_xyz ;;
+    "110") f110_replicate_structure ;;
+    "out2exyz") out2exyz ;;
+    "pos2exyz") pos2extxyz ;;
+    "cif2pos") cif2poscar ;;
+    "cif2exyz") cif2extxyz ;;
+    "xdat2exyz") xdatcar2extxyz ;;
+    "pos2lmp") poscar2lammps ;;
+    "lmp2exyz") lammps2extxyz ;;
+    "traj2exyz") traj2extxyz ;;
     "000") menu; main ;;
 esac
 }
diff --git a/src/f2_sample_structures.sh b/src/f2_sample_structures.sh
index 77a3346..ae52298 100644
--- a/src/f2_sample_structures.sh
+++ b/src/f2_sample_structures.sh
@@ -22,30 +22,14 @@ echo " ---------------------------------------------------"
 
 # Sample structures by pynep
 function f202_pynep_sample_structures(){
-# YELLOW="\033[1;33m"; RESET="\033[0m"
-# echo -en "${YELLOW}"
-echo " +-------------------------------------------------+"
-echo " | PyNEP package is no longer actively maintained  |"
-echo " |       Recommend using function 203 instead      |"
-echo " |  Moreover, qNEP model is not supported for now. |"
 echo " +-------------------------------------------------+"
-# echo -en "${RESET}"
+echo " |  Function 202 (pynep sampling) is deprecated.   |"
+echo " |  PyNEP package is no longer actively maintained.|"
+echo " |  Please use 203) neptrain sampling instead.     |"
+echo " |  If you still insist on pynep, run directly:    |"
+echo " |                 gpumdkit.sh -pynep              |"
 echo " +-------------------------------------------------+"
-echo " |     To use parallel version, please use:        |"
-echo " |            gpumdkit.sh -pynep                   |"
-echo " +-------------------------------------------------+"
-echo " >-------------------------------------------------<"
-echo " | Calling the script in Scripts/sample_structures |"
-echo " | Script: pynep_select_structs.py                 |"
-echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
-echo " >-------------------------------------------------<"
-echo " Input <sample.xyz> <train.xyz> <nep_model>"
-echo " Examp: dump.xyz train.xyz nep.txt "
-echo " ------------>>"
-read -p " " sample_choice
-echo " ---------------------------------------------------"
-python ${GPUMDkit_path}/Scripts/sample_structures/pynep_select_structs.py ${sample_choice}
-echo " Code path: ${GPUMDkit_path}/Scripts/sample_structures/pynep_select_structs.py"
+exit 0
 }
 
 # Parallel sample structures by pynep
@@ -87,9 +71,15 @@ echo " ------------>>"
 read -p " " sample_choice
 echo " ---------------------------------------------------"
 python ${GPUMDkit_path}/Scripts/sample_structures/neptrain_select_structs.py ${sample_choice}
-rm dpdispatcher.log
+rm -f dpdispatcher.log
 echo " Code path: ${GPUMDkit_path}/Scripts/sample_structures/neptrain_select_structs.py"
 echo " ---------------------------------------------------"
+echo " "
+echo " +--------------------------------------------------------+"
+echo " | For this function, please cite the NepTrain paper:     |"
+echo " | Chen et al. Comput. Phys. Commun., 2025, 317, 109859.  |"
+echo " | URL: https://doi.org/10.1016/j.cpc.2025.109859         |"
+echo " +--------------------------------------------------------+"
 }
 
 # Perturb structure
@@ -132,19 +122,22 @@ echo " ---------------------------------------------------"
 
 #--------------------- function 2 ----------------------
 function f2_sample_structures(){
-echo " ------------>>"
-echo " 201) Sample structures from extxyz"
-echo " 202) Sample structures by pynep"
-echo " 203) Sample structures by neptrain"
-echo " 204) Perturb structure"
-echo " 205) Select max force deviation structs"
-echo " 000) Return to the main menu"
-echo " ------------>>"
+echo " +------------------------------------------------------+"
+echo " |                 SAMPLE STRUCTURE TOOLS               |"
+echo " +------------------------------------------------------+"
+echo " | 201) Sample structures from extxyz                   |"
+echo " | 202) FPS sampling by PyNEP [deprecated]              |"
+echo " | 203) FPS sampling by NepTrain [preferred]            |"
+echo " | 204) Perturb structure                               |"
+echo " | 205) Select max force deviation structs              |"
+echo " +------------------------------------------------------+"
+echo " | 000) Return to the main menu                         |"
+echo " +------------------------------------------------------+"
 echo " Input the function number:"
 
-arry_num_choice=("000" "201" "202" "203" "204" "205") 
+valid_menu_choices=("000" "201" "202" "203" "204" "205")
 read -p " " num_choice
-while ! echo "${arry_num_choice[@]}" | grep -wq "$num_choice" 
+while ! echo "${valid_menu_choices[@]}" | grep -wq "$num_choice"
 do
   echo " ------------>>"
   echo " Please reinput function number..."
@@ -159,4 +152,4 @@ case $num_choice in
     "205") f205_select_max_force_deviation_structs ;;
     "000") menu; main ;;
 esac
-}
\ No newline at end of file
+}
diff --git a/src/f3_workflows.sh b/src/f3_workflows.sh
index bfce4b0..dbf9b2d 100644
--- a/src/f3_workflows.sh
+++ b/src/f3_workflows.sh
@@ -34,14 +34,16 @@ echo " ---------------------------------------------------"
 
 # SCF batch pretreatment (main function)
 function f301_scf_batch_pretreatment(){ 
-    echo " ------------>>"
-    echo " 1) VASP scf batch pretreatment"
-    echo " 2) CP2K scf batch pretreatment"
-    echo " ------------>>"
+    echo " +-------------------------------------------------------------+"
+    echo " |                 SCF BATCH PRETREATMENT TOOLS               |"
+    echo " +-------------------------------------------------------------+"
+    echo " | 1) VASP SCF batch pretreatment                              |"
+    echo " | 2) CP2K SCF batch pretreatment                              |"
+    echo " +-------------------------------------------------------------+"
     echo " Input the function number:"
-    arry_num_choice=("1" "2") 
+    valid_menu_choices=("1" "2")
     read -p " " num_choice
-    while ! echo "${arry_num_choice[@]}" | grep -wq "$num_choice" 
+    while ! echo "${valid_menu_choices[@]}" | grep -wq "$num_choice"
     do
       echo " ------------>>"
       echo " Please reinput function number..."
@@ -60,17 +62,20 @@ function f301_scf_batch_pretreatment(){
 
 # main function of workflow
 function f3_workflow_dev(){
-echo " ------------>>"
-echo " 301) SCF batch pretreatment"
-echo " 302) MD sample batch pretreatment (gpumd)"
-echo " 303) MD sample batch pretreatment (lmp)"
-echo " 000) Return to the main menu"
-echo " ------------>>"
+echo " +---------------------------------------------------------+"
+echo " |                      WORKFLOW TOOLS                     |"
+echo " +---------------------------------------------------------+"
+echo " | 301) SCF batch pretreatment                             |"
+echo " | 302) MD sample batch pretreatment (gpumd)               |"
+echo " | 303) MD sample batch pretreatment (lmp)                 |"
+echo " +---------------------------------------------------------+"
+echo " | 000) Return to the main menu                            |"
+echo " +---------------------------------------------------------+"
 echo " Input the function number:"
 
-arry_num_choice=("000" "301" "302" "303") 
+valid_menu_choices=("000" "301" "302" "303")
 read -p " " num_choice
-while ! echo "${arry_num_choice[@]}" | grep -wq "$num_choice" 
+while ! echo "${valid_menu_choices[@]}" | grep -wq "$num_choice"
 do
   echo " ------------>>"
   echo " Please reinput function number..."
@@ -94,4 +99,4 @@ case $num_choice in
         main
         ;;
 esac
-}
\ No newline at end of file
+}
diff --git a/src/f4_calculators.sh b/src/f4_calculators.sh
index 446d040..9980108 100644
--- a/src/f4_calculators.sh
+++ b/src/f4_calculators.sh
@@ -94,6 +94,7 @@ echo " | Script: calc_neighbor_list.py                   |"
 echo " | Developer: Denan LI (lidenan@westlake.edu.cn)   |"
 echo " >-------------------------------------------------<"
 echo " Input script arguments (eg. -c 4 -n 12 -C Pb Sr -E O)"
+echo " See https://gpumdkit.cn/htmls/polar_material_analysis.html for details"
 echo " ------------>>"
 read -p " " input_calc_neighbor_list
 echo " ---------------------------------------------------"
@@ -110,6 +111,7 @@ echo " | Script: calc_displacement.py                    |"
 echo " | Developer: Denan LI (lidenan@westlake.edu.cn)   |"
 echo " >-------------------------------------------------<"
 echo " Input script arguments (eg. -i movie.xyz -n nl-Pb-O.dat -o displacements.dat)"
+echo " See https://gpumdkit.cn/htmls/polar_material_analysis.html for details"
 echo " ------------>>"
 read -p " " input_calc_displacement
 echo " ---------------------------------------------------"
@@ -126,6 +128,7 @@ echo " | Script: calc_averaged_structure.py              |"
 echo " | Developer: Denan LI (lidenan@westlake.edu.cn)   |"
 echo " >-------------------------------------------------<"
 echo " Input script arguments (eg. -i movie.xyz -l 0.2 -o averaged_structure.xyz)"
+echo " See https://gpumdkit.cn/htmls/polar_material_analysis.html for details"
 echo " ------------>>"
 read -p " " input_calc_averaged_structure
 echo " ---------------------------------------------------"
@@ -142,6 +145,7 @@ echo " | Script: calc_oct_tilt.py                        |"
 echo " | Developer: Denan LI (lidenan@westlake.edu.cn)   |"
 echo " >-------------------------------------------------<"
 echo " Input script arguments (eg. -i model.xyz -n nl-Ti-O.dat -o octahedral_tilt.dat)"
+echo " See https://gpumdkit.cn/htmls/polar_material_analysis.html for details"
 echo " ------------>>"
 read -p " " input_calc_oct_tilt
 echo " ---------------------------------------------------"
@@ -158,6 +162,7 @@ echo " | Script: calc_polarization_abo3.py               |"
 echo " | Developer: Denan LI (lidenan@westlake.edu.cn)   |"
 echo " >-------------------------------------------------<"
 echo " Input script arguments (eg. --nl-ba nl-Ti-A.dat --nl-bo nl-Ti-O.dat --bec Pb=2.5 Sr=2.0 Ti=4.0 O=-2.0)"
+echo " See https://gpumdkit.cn/htmls/polar_material_analysis.html for details"
 echo " ------------>>"
 read -p " " input_calc_polarization_abo3
 echo " ---------------------------------------------------"
@@ -166,26 +171,68 @@ echo " Code path: ${GPUMDkit_path}/Scripts/calculators/calc_polarization_abo3.py
 echo " ---------------------------------------------------"
 }
 
-# main function of calculators
-function f4_calculators(){
+function f411_minimize_structure_by_nep(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in calculators   |"
+echo " | Script: minimize_structure_by_nep.py            |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <structure_file> <nep.txt> [fmax=0.01] [max_steps=1000]"
+echo " Supported structure formats: POSCAR, .xyz"
+echo " Optional arguments: "
+echo "  fmax: Force convergence threshold in eV/Ang (default: 0.01)"
+echo "  max_steps: Maximum number of optimization steps (default: 1000)"
+echo " Example: POSCAR nep.txt 0.01 1000"
 echo " ------------>>"
-echo " 401) Calc ionic conductivity"
-echo " 402) Calc properties by nep"
-echo " 403) Calc descriptors of specific elements"
-echo " 404) Calc density of atomistic states (DOAS)"
-echo " 405) Calc nudged elastic band (NEB) by nep"
-echo " 406) Build neighbor list"
-echo " 407) Calc displacement from trajectory"
-echo " 408) Calc averaged structure"
-echo " 409) Calc octahedral tilt"
-echo " 410) Calc polarization for ABO3"
-echo " 000) Return to the main menu"
+read -p " " input_minimize_structure
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/calculators/calc_minimize.py ${input_minimize_structure}
+echo " Code path: ${GPUMDkit_path}/Scripts/calculators/calc_minimize.py"
+echo " ---------------------------------------------------"
+}
+
+function f412_calc_msd_from_trajectory(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in calculators   |"
+echo " | Script: calc_msd.py                             |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <extxyz_file> <element_symbol> <dt_fs> [max_corr_steps]"
+echo "   Optional argument: max_corr_steps (default: frame number)"
+echo " Example: dump.xyz Li 10"
 echo " ------------>>"
+read -p " " input_calc_msd
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/calculators/calc_msd.py ${input_calc_msd}
+echo " Code path: ${GPUMDkit_path}/Scripts/calculators/calc_msd.py"
+echo " ---------------------------------------------------"
+}
+
+# main function of calculators
+function f4_calculators(){
+echo " +----------------------------------------------------------+"
+echo " |                     CALCULATOR TOOLS                     |"
+echo " +----------------------------------------------------------+"
+echo " | 401) Calc ionic conductivity                             |"
+echo " | 402) Calc properties by nep                              |"
+echo " | 403) Calc descriptors of specific elements               |"
+echo " | 404) Calc density of atomistic states (DOAS)             |"
+echo " | 405) Calc nudged elastic band (NEB) by nep               |"
+echo " | 406) Build neighbor list                                 |"
+echo " | 407) Calc displacement from trajectory                   |"
+echo " | 408) Calc averaged structure                             |"
+echo " | 409) Calc octahedral tilt                                |"
+echo " | 410) Calc polarization for ABO3                          |"
+echo " | 411) Minimize structure by nep                           |"
+echo " | 412) Calc mean square displacement (MSD) from trajectory |"
+echo " +----------------------------------------------------------+"
+echo " | 000) Return to the main menu                             |"
+echo " +----------------------------------------------------------+"
 echo " Input the function number:"
 
-arry_num_choice=("000" "401" "402" "403" "404" "405" "406" "407" "408" "409" "410")
+valid_menu_choices=("000" "401" "402" "403" "404" "405" "406" "407" "408" "409" "410" "411" "412")
 read -p " " num_choice
-while ! echo "${arry_num_choice[@]}" | grep -wq "$num_choice" 
+while ! echo "${valid_menu_choices[@]}" | grep -wq "$num_choice"
 do
   echo " ------------>>"
   echo " Please reinput function number..."
@@ -203,6 +250,8 @@ case $num_choice in
     "408") f408_calc_averaged_structure ;;
     "409") f409_calc_oct_tilt ;;
     "410") f410_calc_polarization_abo3 ;;
+    "411") f411_minimize_structure_by_nep ;;
+    "412") f412_calc_msd_from_trajectory ;;
     "000") menu; main ;;
 esac
 }
diff --git a/src/f5_analyzers.sh b/src/f5_analyzers.sh
index b16c648..5d10ad0 100644
--- a/src/f5_analyzers.sh
+++ b/src/f5_analyzers.sh
@@ -10,7 +10,7 @@ echo " | Script: analyze_composition.py                  |"
 echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
 echo " >-------------------------------------------------<"
 echo " Input <input.xyz> you want to analyze"
-echo " Examp: train.xyz"
+echo " Example: train.xyz"
 echo " ------------>>"
 read -p " " input_extxyz
 echo " ---------------------------------------------------"
@@ -33,16 +33,142 @@ echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/find_outliers.py"
 echo " ---------------------------------------------------"
 }
 
-# main function of analyzers
-function f5_analyzers(){
+# Analyze chemical species of extxyz
+function f503_analyze_chem_species(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in analyzer      |"
+echo " | Script: analyze_chem_species.py                 |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> you want to analyze"
+echo " Example: train.xyz"
+echo " ------------>>"
+read -p " " input_extxyz
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/analyzer/analyze_chem_species.py ${input_extxyz}
+echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/analyze_chem_species.py"
+echo " ---------------------------------------------------"
+}
+
+# Check charge balance of extxyz
+function f504_charge_balance_check(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in analyzer      |"
+echo " | Script: charge_balance_check.py                 |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> you want to analyze"
+echo " Examp: train.xyz"
 echo " ------------>>"
-echo " 501) Analyze composition of extxyz"
-echo " 502) Find outliers of extxyz"
-echo " 000) Return to the main menu"
+read -p " " input_extxyz
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/analyzer/charge_balance_check.py ${input_extxyz}
+echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/charge_balance_check.py"
+echo " ---------------------------------------------------"
+}
+
+# Analyze energy/force/virial range
+function f505_energy_force_virial_analyzer(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in analyzer      |"
+echo " | Script: energy_force_virial_analyzer.py         |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> and property (energy/force/virial)"
+echo " Example: train.xyz force"
+echo " ------------>>"
+read -p " " input_extxyz property_name
+echo " Plot histogram? (y/n)"
 echo " ------------>>"
+read -p " " if_hist
+echo " ---------------------------------------------------"
+if [ "$if_hist" == "y" ] || [ "$if_hist" == "Y" ]; then
+  python ${GPUMDkit_path}/Scripts/analyzer/energy_force_virial_analyzer.py ${input_extxyz} ${property_name} hist
+else
+  python ${GPUMDkit_path}/Scripts/analyzer/energy_force_virial_analyzer.py ${input_extxyz} ${property_name}
+fi
+echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/energy_force_virial_analyzer.py"
+echo " ---------------------------------------------------"
+}
+
+# Filter structures by minimum distance (without PBC)
+function f506_filter_structures_by_distance(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in analyzer      |"
+echo " | Script: filter_structures_by_distance.py        |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> and <min_dist>"
+echo " Example: dump.xyz 1.5"
+echo " ------------>>"
+read -p " " input_extxyz min_dist
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/analyzer/filter_structures_by_distance.py ${input_extxyz} ${min_dist}
+echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/filter_structures_by_distance.py"
+echo " ---------------------------------------------------"
+}
+
+# Get minimum interatomic distance with optional PBC
+function f507_get_min_dist(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in analyzer      |"
+echo " | Script: get_min_dist.py / get_min_dist_pbc.py   |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <input.xyz> you want to analyze"
+echo " Example: dump.xyz"
+echo " ------------>>"
+read -p " " input_extxyz
+echo " Consider PBC? (y/n)"
+echo " ------------>>"
+read -p " " consider_pbc
+echo " ---------------------------------------------------"
+if [ "$consider_pbc" == "y" ] || [ "$consider_pbc" == "Y" ]; then
+  python ${GPUMDkit_path}/Scripts/analyzer/get_min_dist_pbc.py ${input_extxyz}
+  echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/get_min_dist_pbc.py"
+else
+  python ${GPUMDkit_path}/Scripts/analyzer/get_min_dist.py ${input_extxyz}
+  echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/get_min_dist.py"
+fi
+echo " ---------------------------------------------------"
+}
+
+# Probability density analysis for diffusion channels
+function f508_probability_density_analysis(){
+echo " >-------------------------------------------------<"
+echo " | This function calls the script in analyzer      |"
+echo " | Script: probability_density_analysis.py         |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input <ref_struct> <trajectory_file> <species> <interval>"
+echo " Example: LLZO.vasp dump.xyz Li 0.25"
+echo " ------------>>"
+read -p " " ref_struct trajectory_file species interval
+echo " ---------------------------------------------------"
+python ${GPUMDkit_path}/Scripts/analyzer/probability_density_analysis.py ${ref_struct} ${trajectory_file} ${species} ${interval}
+echo " Code path: ${GPUMDkit_path}/Scripts/analyzer/probability_density_analysis.py"
+echo " ---------------------------------------------------"
+}
+
+# main function of analyzers
+function f5_analyzers(){
+echo " +------------------------------------------------------+"
+echo " |                    ANALYZER TOOLS                    |"
+echo " +------------------------------------------------------+"
+echo " | 501) Analyze composition of extxyz                   |"
+echo " | 502) Find outliers of extxyz                         |"
+echo " | 503) Analyze chemical species of extxyz              |"
+echo " | 504) Check charge balance of extxyz                  |"
+echo " | 505) Analyze energy/force/virial range               |"
+echo " | 506) Filter structures by minimum distance           |"
+echo " | 507) Get minimum interatomic distance                |"
+echo " | 508) Probability density analysis                    |"
+echo " +------------------------------------------------------+"
+echo " | 000) Return to the main menu                         |"
+echo " +------------------------------------------------------+"
 echo " Input the function number:"
 
-arry_num_choice=("000" "501" "502") 
+arry_num_choice=("000" "501" "502" "503" "504" "505" "506" "507" "508") 
 read -p " " num_choice
 while ! echo "${arry_num_choice[@]}" | grep -wq "$num_choice" 
 do
@@ -54,6 +180,12 @@ done
 case $num_choice in
     "501") f501_analyze_composition ;;            
     "502") f502_find_outliers ;;
+    "503") f503_analyze_chem_species ;;
+    "504") f504_charge_balance_check ;;
+    "505") f505_energy_force_virial_analyzer ;;
+    "506") f506_filter_structures_by_distance ;;
+    "507") f507_get_min_dist ;;
+    "508") f508_probability_density_analysis ;;
     "000") menu; main ;;
 esac
-}
\ No newline at end of file
+}
diff --git a/src/f6_plots.sh b/src/f6_plots.sh
new file mode 100644
index 0000000..27eb47e
--- /dev/null
+++ b/src/f6_plots.sh
@@ -0,0 +1,107 @@
+function f6_plots_one_column(){
+    echo " +---------------------------------------------------------+"
+    echo " |                PLOT & VISUALIZATION TOOLS               |"
+    echo " +---------------------------------------------------------+"
+    echo " |  Usage: gpumdkit.sh -plt <type>                         |"
+    echo " +---------------------------------------------------------+"
+    echo " |                NEP Training & Evaluation                |"
+    echo " +---------------------------------------------------------+"
+    echo " |  train          Plot NEP training results               |"
+    echo " |  prediction     Plot NEP prediction results             |"
+    echo " |  train_test     Plot NEP train and test results         |"
+    echo " |  parity_density Plot parity density plot                |"
+    echo " |  train_density  Plot training results in density plot   |" 
+    echo " |  restart        Plot parameters in nep.restart          |"
+    echo " |  charge         Plot charge distribution                |"
+    echo " |  born_charge    Plot Born effective charges             |"    
+    echo " |  dimer          Plot dimer energy/force curve           |"
+    echo " |  force_errors   Plot force errors                       |"
+    echo " |  des            Plot descriptors                        |"       
+    echo " |  lr             Plot learning rate for gnep             |"
+    echo " +---------------------------------------------------------+"
+    echo " |                 Diffusion & Transport                   |"
+    echo " +---------------------------------------------------------+"
+    echo " |  msd            Plot mean square displacement (MSD)     |"
+    echo " |  msd_conv       Plot MSD convergence                    |"
+    echo " |  msd_all        Plot MSD for all species                |"
+    echo " |  sdc            Plot self diffusion coefficient         |"
+    echo " |  msd_sdc        Plot MSD and SDC together               |"
+    echo " |  sigma          Plot Arrhenius ionic conductivity       |"
+    echo " |  D              Plot Arrhenius diffusivity              |"
+    echo " |  sigma_xyz      Plot directional Arrhenius sigma        |"
+    echo " |  D_xyz          Plot directional Arrhenius D            |"
+    echo " |  doas           Plot density of atomistic states        |"
+    echo " +---------------------------------------------------------+"
+    echo " |                 MD & Structural Analysis                |"
+    echo " +---------------------------------------------------------+"
+    echo " |  thermo         Plot thermo info in thermo.out          |"
+    echo " |  thermo2/3      Plot thermo in different styles         |"    
+    echo " |  rdf            Plot radial distribution function       |"
+    echo " |  rdf_pmf        Plot potential of mean force            |"
+    echo " |  vac            Plot velocity autocorrelation           |"
+    echo " |  cohesive       Plot cohesive energy curve              |"
+    echo " |  net_force      Plot net force distribution             |"
+    echo " |  plane-grid     Plot displacement plane grid            |"
+    echo " +---------------------------------------------------------+"
+    echo " |                     Heat Transport                      |"
+    echo " +---------------------------------------------------------+"
+    echo " |  emd            Plot EMD results                        |"
+    echo " |  nemd           Plot NEMD results                       |"
+    echo " |  hnemd          Plot HNEMD results                      |"
+    echo " |  viscosity      Plot viscosity                          |"
+    echo " +---------------------------------------------------------+"
+    echo " |                         Phonons                         |"
+    echo " +---------------------------------------------------------+"
+    echo " |  pdos           Plot VAC and PDOS                       |"
+    echo " +---------------------------------------------------------+"
+    # echo " |                   Miscellaneous Tools                   |"
+    # echo " +---------------------------------------------------------+"
+    # echo " |  Developing ...                                         |"
+    # echo " +---------------------------------------------------------+"
+}
+
+function f6_plots_two_column(){
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |                     GPUMDkit ${VERSION} PLOT & VISUALIZATION TOOLS              |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |  Usage: gpumdkit.sh -plt <type>                        Help: gpumdkit.sh -plt <type> -h       |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |                                    NEP Training & Evaluation                                  |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |  train          - NEP training results           prediction     - NEP prediction results      |"
+    echo " |  train_test     - NEP train and test results     parity_density - Parity density plot         |"
+    echo " |  train_density  - Training results density plot  restart        - Parameters in nep.restart   |"
+    echo " |  charge         - Charge distribution            born_charge    - Born effective charges      |"
+    echo " |  dimer          - Dimer energy/force curve       force_errors   - Force errors                |"
+    echo " |  des            - Descriptors                    lr             - Learning rate for gnep      |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |                                     Diffusion & Transport                                     |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |  msd            - Mean square displacement       msd_conv       - MSD convergence             |"
+    echo " |  msd_all        - MSD for all species            sdc            - Self diffusion coefficient  |"
+    echo " |  msd_sdc        - MSD and SDC together           sigma          - Arrhenius ionic conductivity|"
+    echo " |  D              - Arrhenius diffusivity          sigma_xyz      - Directional Arrhenius sigma |"
+    echo " |  D_xyz          - Directional Arrhenius D                                                     |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |                                    MD & Structural Analysis                                   |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |  thermo         - thermo info in thermo.out      thermo2/3      - Thermo in different styles  |"
+    echo " |  rdf            - Radial distribution function   rdf_pmf        - Potential of mean force     |"
+    echo " |  vac            - Velocity autocorrelation       cohesive       - Cohesive energy curve       |"
+    echo " |  net_force      - Net force distribution         plane-grid     - Displacement plane grid     |"
+    echo " |  doas           - Density of atomistic states                                                 |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |                                        Heat Transport                                         |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |  emd            - EMD results                    nemd           - NEMD results                |"
+    echo " |  hnemd          - HNEMD results                  viscosity      - Viscosity                   |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |                                          Phonons                                              |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    echo " |  pdos           - VAC and PDOS                                                                |"
+    echo " +-----------------------------------------------------------------------------------------------+"
+    # echo " |                                       Miscellaneous Tools                                     |"
+    # echo " +-----------------------------------------------------------------------------------------------+"
+    # echo " |  Developing ...                                                                               |"
+    # echo " +-----------------------------------------------------------------------------------------------+"
+}
diff --git a/src/f7_utilities.sh b/src/f7_utilities.sh
new file mode 100644
index 0000000..81f6130
--- /dev/null
+++ b/src/f7_utilities.sh
@@ -0,0 +1,48 @@
+#--------------------- function 7 utilities ----------------------
+# Miscellaneous tools that do not fit other modules
+
+function f701_time_consuming_analyzer(){
+echo " >-------------------------------------------------<"
+echo " | Calling scripts in Scripts/analyzer             |"
+echo " | Script: time_consuming_gpumd.sh/nep.sh          |"
+echo " | Developer: Zihan YAN (yanzihan@westlake.edu.cn) |"
+echo " >-------------------------------------------------<"
+echo " Input mode: gpumd, nep or gnep"
+echo " Example: gpumd"
+echo " ------------>>"
+read -p " " time_mode
+echo " ---------------------------------------------------"
+case $time_mode in
+    gpumd) bash ${analyzer_path}/time_consuming_gpumd.sh ;;
+    nep|gnep) bash ${analyzer_path}/time_consuming_nep.sh ;;
+    *)
+        echo " Please input 'gpumd', 'nep' or 'gnep'." ;;
+esac
+echo " Code path: ${analyzer_path}/time_consuming_*.sh"
+echo " ---------------------------------------------------"
+}
+
+function f7_utilities(){
+echo " +------------------------------------------------------+"
+echo " |                    UTILITY TOOLS                     |"
+echo " +------------------------------------------------------+"
+echo " | 701) Time consuming analyzer (gpumd/nep)             |"
+echo " +------------------------------------------------------+"
+echo " | 000) Return to the main menu                         |"
+echo " +------------------------------------------------------+"
+echo " Input the function number:"
+
+valid_menu_choices=("000" "701")
+read -p " " num_choice
+while ! echo "${valid_menu_choices[@]}" | grep -wq "$num_choice"
+do
+  echo " ------------>>"
+  echo " Please reinput function number..."
+  read -p " " num_choice
+done
+
+case $num_choice in
+    "701") f701_time_consuming_analyzer ;;
+    "000") menu; main ;;
+esac
+}