diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index eef096d..f0fd20b 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -3,7 +3,6 @@ repos:
     rev: v4.2.0
     hooks:
       - id: end-of-file-fixer
-      - id: fix-encoding-pragma
       - id: mixed-line-ending
       - id: trailing-whitespace
       - id: check-json
@@ -24,4 +23,4 @@ repos:
     rev: 6.0.0
     hooks:
       - id: flake8
-        args: ["--max-line-length=88", "--ignore=E203,W503"]
\ No newline at end of file
+        args: ["--max-line-length=88", "--ignore=E203,W503"]
diff --git a/src/mattersim/applications/relax.py b/src/mattersim/applications/relax.py
index f18ed27..9a678ac 100644
--- a/src/mattersim/applications/relax.py
+++ b/src/mattersim/applications/relax.py
@@ -8,6 +8,7 @@
 from ase.optimize import BFGS, FIRE
 from ase.optimize.optimize import Optimizer
 from ase.units import GPa
+from deprecated import deprecated
 
 
 class Relaxer(object):
@@ -55,7 +56,7 @@ def relax(
         fmax: float = 0.01,
         params_filter: dict = {},
         **kwargs,
-    ) -> Atoms:
+    ) -> Tuple[bool, Atoms]:
         """
         Relax the atoms object.
 
@@ -108,6 +109,7 @@ def relax(
         return converged, atoms
 
     @classmethod
+    @deprecated(reason="Use cli/applications/relax_structure.py instead.")
     def relax_structures(
         cls,
         atoms: Union[Atoms, Iterable[Atoms]],
diff --git a/src/mattersim/cli/__init__.py b/src/mattersim/cli/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/mattersim/cli/applications/__init__.py b/src/mattersim/cli/applications/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/mattersim/cli/applications/moldyn.py b/src/mattersim/cli/applications/moldyn.py
new file mode 100644
index 0000000..a6b8219
--- /dev/null
+++ b/src/mattersim/cli/applications/moldyn.py
@@ -0,0 +1,114 @@
+import os
+import re
+import uuid
+from collections import defaultdict
+from typing import List
+
+import pandas as pd
+from ase import Atoms
+from ase.io import read
+from loguru import logger
+from pymatgen.io.ase import AseAtomsAdaptor
+
+from mattersim.applications.moldyn import MolecularDynamics
+
+
+def moldyn(
+    atoms_list: List[Atoms],
+    *,
+    temperature: float = 300,
+    timestep: float = 1,
+    steps: int = 1000,
+    ensemble: str = "nvt_nose_hoover",
+    logfile: str = "-",
+    loginterval: int = 10,
+    trajectory: str = None,
+    taut: float = None,
+    work_dir: str = str(uuid.uuid4()),
+    save_csv: str = "results.csv.gz",
+    **kwargs,
+) -> dict:
+    moldyn_results = defaultdict(list)
+
+    for atoms in atoms_list:
+        # check if the atoms object has non-zero values in the lower triangle
+        # of the cell. If so, the cell will be rotated and permuted to upper
+        # triangular form. This is to avoid numerical issues in the MD simulation.
+        print(atoms.cell.array)
+        if any(atoms.cell.array[2, 0:2]) or atoms.cell.array[1, 0] != 0:
+            logger.warning(
+                "The lower triangle of the cell is not zero. "
+                "The cell will be rotated and permuted to upper triangular form."
+            )
+
+            # The following code is from the PR
+            # https://gitlab.com/ase/ase/-/merge_requests/3277.
+            # It will be removed once the PR is merged.
+            # This part of the codes rotates the cell and permutes the axes
+            # such that the cell will be in upper triangular form.
+
+            from ase.build import make_supercell
+
+            _calc = atoms.calc
+            logger.info(f"Initial cell: {atoms.cell.array}")
+
+            atoms.set_cell(atoms.cell.standard_form()[0], scale_atoms=True)
+
+            # Permute a and c axes
+            atoms = make_supercell(atoms, [[0, 0, 1], [0, 1, 0], [1, 0, 0]])
+
+            atoms.rotate(90, "y", rotate_cell=True)
+
+            # set the lower triangle of the cell to be exactly zero
+            # to avoid numerical issues
+            atoms.cell.array[1, 0] = 0
+            atoms.cell.array[2, 0] = 0
+            atoms.cell.array[2, 1] = 0
+
+            logger.info(f"Cell after rotation/permutation: {atoms.cell.array}")
+            atoms.calc = _calc
+
+        if not os.path.exists(work_dir):
+            os.makedirs(work_dir)
+
+        md = MolecularDynamics(
+            atoms,
+            ensemble=ensemble,
+            temperature=temperature,
+            timestep=timestep,
+            loginterval=loginterval,
+            logfile=os.path.join(work_dir, logfile),
+            trajectory=os.path.join(work_dir, trajectory),
+            taut=taut,
+        )
+        md.run(steps)
+
+        # parse the logfile
+
+        # Read the file into a pandas DataFrame
+        df = pd.read_csv(
+            os.path.join(work_dir, logfile),
+            sep="\\s+",
+            names=["time", "temperature", "energy", "pressure"],
+            skipfooter=1,
+        )
+        df.columns = list(
+            map(lambda x: re.sub(r"\[.*?\]", "", x).strip().lower(), df.columns)
+        )
+        traj = read(os.path.join(work_dir, trajectory), index=":")
+        print(df.shape)
+        print(len(traj))
+        structure_list = [AseAtomsAdaptor.get_structure(atoms) for atoms in traj]
+
+        # Add the structure list to the DataFrame
+        df["structure"] = [structure.to_json() for structure in structure_list]
+
+        # Print the DataFrame
+        print(df)
+
+        # Save the DataFrame to a CSV file
+        df.to_csv(os.path.join(work_dir, save_csv))
+
+        moldyn_results = df.to_dict()
+
+    return moldyn_results
diff --git a/src/mattersim/cli/applications/phonon.py b/src/mattersim/cli/applications/phonon.py
new file mode 100644
index 0000000..090cc81
--- /dev/null
+++ b/src/mattersim/cli/applications/phonon.py
@@ -0,0 +1,142 @@
+import os
+import uuid
+from collections import defaultdict
+from typing import List
+
+import numpy as np
+import pandas as pd
+import yaml
+from ase import Atoms
+from loguru import logger
+from pymatgen.core.structure import Structure
+from pymatgen.io.ase import AseAtomsAdaptor
+from tqdm import tqdm
+
+from mattersim.applications.phonon import PhononWorkflow
+from mattersim.cli.applications.relax import relax
+
+
+def phonon(
+    atoms_list: List[Atoms],
+    *,
+    find_prim: bool = False,
+    work_dir: str = str(uuid.uuid4()),
+    save_csv: str = "results.csv.gz",
+    amplitude: float = 0.01,
+    supercell_matrix: np.ndarray = None,
+    qpoints_mesh: np.ndarray = None,
+    max_atoms: int = None,
+    enable_relax: bool = False,
+    **kwargs,
+) -> dict:
+    """
+    Predict phonon properties for a list of atoms.
+
+    Args:
+        atoms_list (List[Atoms]): List of ASE Atoms objects.
+        find_prim (bool, optional): If find the primitive cell and use it
+            to calculate phonon. Default to False.
+        work_dir (str, optional): workplace path to contain phonon result.
+            Defaults to data + chemical_symbols + 'phonon'
+        amplitude (float, optional): Magnitude of the finite difference to
+            displace in force constant calculation, in Angstrom. Defaults
+            to 0.01 Angstrom.
+        supercell_matrix (nd.array, optional): Supercell matrix for constr
+            -uct supercell, priority over than max_atoms. Defaults to None.
+        qpoints_mesh (nd.array, optional): Qpoint mesh for IBZ integral,
+            priority over than max_atoms. Defaults to None.
+        max_atoms (int, optional): Maximum atoms number limitation for the
+            supercell generation. If not set, will automatic generate super
+            -cell based on symmetry. Defaults to None.
+        enable_relax (bool, optional): Whether to relax the structure before
+            predicting phonon properties. Defaults to False.
+    """
+    phonon_results = defaultdict(list)
+
+    for atoms in tqdm(
+        atoms_list, total=len(atoms_list), desc="Predicting phonon properties"
+    ):
+        if enable_relax:
+            relaxed_results = relax(
+                [atoms],
+                constrain_symmetry=True,
+                work_dir=work_dir,
+                save_csv=save_csv.replace(".csv", "_relax.csv"),
+            )
+            structure = Structure.from_str(relaxed_results["structure"][0], fmt="json")
+            _atoms = AseAtomsAdaptor.get_atoms(structure)
+            _atoms.calc = atoms.calc
+            atoms = _atoms
+        ph = PhononWorkflow(
+            atoms=atoms,
+            find_prim=find_prim,
+            work_dir=work_dir,
+            amplitude=amplitude,
+            supercell_matrix=supercell_matrix,
+            qpoints_mesh=qpoints_mesh,
+            max_atoms=max_atoms,
+        )
+        has_imaginary, phonon = ph.run()
+        phonon_results["has_imaginary"].append(has_imaginary)
+        # phonon_results["phonon"].append(phonon)
+        phonon_results["phonon_band_plot"].append(
+            os.path.join(os.path.abspath(work_dir), f"{atoms.symbols}_phonon_band.png")
+        )
+        phonon_results["phonon_dos_plot"].append(
+            os.path.join(os.path.abspath(work_dir), f"{atoms.symbols}_phonon_dos.png")
+        )
+        os.rename(
+            os.path.join(os.path.abspath(work_dir), "band.yaml"),
+            os.path.join(os.path.abspath(work_dir), f"{atoms.symbols}_band.yaml"),
+        )
+        os.rename(
+            os.path.join(os.path.abspath(work_dir), "phonopy_params.yaml"),
+            os.path.join(
+                os.path.abspath(work_dir), f"{atoms.symbols}_phonopy_params.yaml"
+            ),
+        )
+        os.rename(
+            os.path.join(os.path.abspath(work_dir), "total_dos.dat"),
+            os.path.join(os.path.abspath(work_dir), f"{atoms.symbols}_total_dos.dat"),
+        )
+        phonon_results["phonon_band"].append(
+            yaml.safe_load(
+                open(
+                    os.path.join(
+                        os.path.abspath(work_dir), f"{atoms.symbols}_band.yaml"
+                    ),
+                    "r",
+                )
+            )
+        )
+        phonon_results["phonopy_params"].append(
+            yaml.safe_load(
+                open(
+                    os.path.join(
+                        os.path.abspath(work_dir),
+                        f"{atoms.symbols}_phonopy_params.yaml",
+                    ),
+                    "r",
+                )
+            )
+        )
+        phonon_results["total_dos"].append(
+            np.loadtxt(
+                os.path.join(
+                    os.path.abspath(work_dir), f"{atoms.symbols}_total_dos.dat"
+                ),
+                comments="#",
+            )
+        )
+
+    if not os.path.exists(work_dir):
+        os.makedirs(work_dir)
+
+    logger.info(f"Saving the results to {os.path.join(work_dir, save_csv)}")
+    df = pd.DataFrame(phonon_results)
+    df.to_csv(
+        os.path.join(work_dir, save_csv.replace(".csv", "_phonon.csv")),
+        index=False,
+        mode="a",
+    )
+    return phonon_results
diff --git a/src/mattersim/cli/applications/relax.py b/src/mattersim/cli/applications/relax.py
new file mode 100644
index 0000000..d477f45
--- /dev/null
+++ b/src/mattersim/cli/applications/relax.py
@@ -0,0 +1,98 @@
+import os
+import uuid
+from collections import defaultdict
+from typing import List, Union
+
+import pandas as pd
+from ase import Atoms
+from ase.constraints import Filter
+from ase.optimize.optimize import Optimizer
+from ase.units import GPa
+from loguru import logger
+from pymatgen.io.ase import AseAtomsAdaptor
+from tqdm import tqdm
+
+from mattersim.applications.relax import Relaxer
+
+
+def relax(
+    atoms_list: List[Atoms],
+    *,
+    optimizer: Union[str, Optimizer] = "FIRE",
+    filter: Union[str, Filter, None] = None,
+    constrain_symmetry: bool = False,
+    fix_axis: Union[bool, List[bool]] = False,
+    pressure_in_GPa: float = None,
+    fmax: float = 0.01,
+    steps: int = 500,
+    work_dir: str = str(uuid.uuid4()),
+    save_csv: str = "results.csv.gz",
+    **kwargs,
+) -> dict:
+    """
+    Relax a list of atoms structures.
+
+    Args:
+        atoms_list (List[Atoms]): List of ASE Atoms objects.
+        optimizer (Union[str, Optimizer]): The optimizer to use. Default is "FIRE".
+        filter (Union[str, Filter, None]): The filter to use.
+        constrain_symmetry (bool): Whether to constrain symmetry. Default is False.
+        fix_axis (Union[bool, List[bool]]): Whether to fix the axis. Default is False.
+        pressure_in_GPa (float): Pressure in GPa to use for relaxation.
+        fmax (float): Maximum force tolerance for relaxation. Default is 0.01.
+        steps (int): Maximum number of steps for relaxation. Default is 500.
+        work_dir (str): Working directory for the calculations.
+            Default is a UUID with timestamp.
+        save_csv (str): Save the results to a CSV file. Default is `results.csv.gz`.
+
+    Returns:
+        pd.DataFrame: DataFrame containing the relaxed results.
+    """
+    params_filter = {}
+
+    if pressure_in_GPa:
+        params_filter["scalar_pressure"] = (
+            pressure_in_GPa * GPa
+        )  # convert GPa to eV/Angstrom^3
+        filter = "ExpCellFilter" if filter is None else filter
+    elif filter:
+        params_filter["scalar_pressure"] = 0.0
+
+    relaxer = Relaxer(
+        optimizer=optimizer,
+        filter=filter,
+        constrain_symmetry=constrain_symmetry,
+        fix_axis=fix_axis,
+    )
+
+    relaxed_results = defaultdict(list)
+    for atoms in tqdm(atoms_list, total=len(atoms_list), desc="Relaxing structures"):
+        converged, relaxed_atoms = relaxer.relax(
+            atoms,
+            params_filter=params_filter,
+            fmax=fmax,
+            steps=steps,
+        )
+        relaxed_results["converged"].append(converged)
+        relaxed_results["structure"].append(
+            AseAtomsAdaptor.get_structure(relaxed_atoms).to_json()
+        )
+        relaxed_results["energy"].append(relaxed_atoms.get_potential_energy())
+        relaxed_results["energy_per_atom"].append(
+            relaxed_atoms.get_potential_energy() / len(relaxed_atoms)
+        )
+        relaxed_results["forces"].append(relaxed_atoms.get_forces())
+        relaxed_results["stress"].append(relaxed_atoms.get_stress(voigt=False))
+        relaxed_results["stress_GPa"].append(
+            relaxed_atoms.get_stress(voigt=False) / GPa
+        )
+
+        logger.info(f"Relaxed structure: {relaxed_atoms}")
+
+    if not os.path.exists(work_dir):
+        os.makedirs(work_dir)
+
+    logger.info(f"Saving the results to {os.path.join(work_dir, save_csv)}")
+    df = pd.DataFrame(relaxed_results)
+    df.to_csv(os.path.join(work_dir, save_csv), index=False, mode="a")
+    return relaxed_results
diff --git a/src/mattersim/cli/applications/singlepoint.py b/src/mattersim/cli/applications/singlepoint.py
new file mode 100644
index 0000000..9c0153f
--- /dev/null
+++ b/src/mattersim/cli/applications/singlepoint.py
@@ -0,0 +1,46 @@
+import os
+import uuid
+from collections import defaultdict
+from typing import List
+
+import pandas as pd
+from ase import Atoms
+from ase.units import GPa
+from loguru import logger
+from pymatgen.io.ase import AseAtomsAdaptor
+from tqdm import tqdm
+
+
+def singlepoint(
+    atoms_list: List[Atoms],
+    *,
+    work_dir: str = str(uuid.uuid4()),
+    save_csv: str = "results.csv.gz",
+    **kwargs,
+) -> dict:
+    """
+    Predict single point properties for a list of atoms.
+
+    """
+    logger.info("Predicting single point properties.")
+    predicted_properties = defaultdict(list)
+    for atoms in tqdm(
+        atoms_list, total=len(atoms_list), desc="Predicting single point properties"
+    ):
+        predicted_properties["structure"].append(AseAtomsAdaptor.get_structure(atoms))
+        predicted_properties["energy"].append(atoms.get_potential_energy())
+        predicted_properties["energy_per_atom"].append(
+            atoms.get_potential_energy() / len(atoms)
+        )
+        predicted_properties["forces"].append(atoms.get_forces())
+        predicted_properties["stress"].append(atoms.get_stress(voigt=False))
+        predicted_properties["stress_GPa"].append(atoms.get_stress(voigt=False) / GPa)
+
+    if not os.path.exists(work_dir):
+        os.makedirs(work_dir)
+
+    logger.info(f"Saving the results to {os.path.join(work_dir, save_csv)}")
+
+    df = pd.DataFrame(predicted_properties)
+    df.to_csv(os.path.join(work_dir, save_csv), index=False, mode="a")
+    return predicted_properties
diff --git a/src/mattersim/cli/mattersim_app.py b/src/mattersim/cli/mattersim_app.py
new file mode 100644
index 0000000..ae51e04
--- /dev/null
+++ b/src/mattersim/cli/mattersim_app.py
@@ -0,0 +1,344 @@
+import argparse
+import uuid
+from datetime import datetime
+from typing import List, Union
+
+from ase import Atoms
+from ase.io import read as ase_read
+from loguru import logger
+
+from mattersim.cli.applications.moldyn import moldyn
+from mattersim.cli.applications.phonon import phonon
+from mattersim.cli.applications.relax import relax
+from mattersim.cli.applications.singlepoint import singlepoint
+from mattersim.forcefield import MatterSimCalculator
+
+
+def singlepoint_cli(args: argparse.Namespace) -> dict:
+    """
+    CLI wrapper for singlepoint function.
+
+    Args:
+        args (argparse.Namespace): Command line arguments.
+
+    Returns:
+        dict: Dictionary containing the predicted properties.
+
+    """
+    atoms_list = parse_atoms_list(
+        args.structure_file, args.mattersim_model, args.device
+    )
+    singlepoint_args = {
+        k: v
+        for k, v in vars(args).items()
+        if k not in ["structure_file", "mattersim_model", "device"]
+    }
+    return singlepoint(atoms_list, **singlepoint_args)
+
+
+def relax_cli(args: argparse.Namespace) -> dict:
+    """
+    CLI wrapper for relax function.
+
+    Args:
+        args (argparse.Namespace): Command line arguments.
+
+    Returns:
+        dict: Dictionary containing the relaxed results
+    """
+    atoms_list = parse_atoms_list(
+        args.structure_file, args.mattersim_model, args.device
+    )
+    relax_args = {
+        k: v
+        for k, v in vars(args).items()
+        if k not in ["structure_file", "mattersim_model", "device"]
+    }
+    return relax(atoms_list, **relax_args)
+
+
+def phonon_cli(args: argparse.Namespace) -> dict:
+    """
+    CLI wrapper for phonon function.
+
+    Args:
+        args (argparse.Namespace): Command line arguments.
+
+    Returns:
+        dict: Dictionary containing the phonon properties.
+    """
+    atoms_list = parse_atoms_list(
+        args.structure_file, args.mattersim_model, args.device
+    )
+    phonon_args = {
+        k: v
+        for k, v in vars(args).items()
+        if k not in ["structure_file", "mattersim_model", "device"]
+    }
+    return phonon(atoms_list, **phonon_args)
+
+
+def moldyn_cli(args: argparse.Namespace) -> dict:
+    """
+    CLI wrapper for moldyn function.
+
+    Args:
+        args (argparse.Namespace): Command line arguments.
+
+    Returns:
+        dict: Dictionary containing the molecular dynamics properties.
+    """
+    atoms_list = parse_atoms_list(
+        args.structure_file, args.mattersim_model, args.device
+    )
+    if len(atoms_list) > 1:
+        logger.error("Molecular dynamics may take too long for multiple structures.")
+
+    moldyn_args = {
+        k: v
+        for k, v in vars(args).items()
+        if k not in ["structure_file", "mattersim_model", "device"]
+    }
+    return moldyn(atoms_list, **moldyn_args)
+
+
+def add_common_args(parser: argparse.ArgumentParser):
+    parser.add_argument(
+        "--structure-file",
+        type=str,
+        nargs="+",
+        help="Path to the atoms structure file(s).",
+    )
+    parser.add_argument(
+        "--mattersim-model",
+        type=str,
+        choices=["mattersim-v1.0.0-1m", "mattersim-v1.0.0-5m"],
+        default="mattersim-v1.0.0-1m",
+        help="MatterSim model to use.",
+    )
+    parser.add_argument(
+        "--device",
+        type=str,
+        default="cpu",
+        choices=["cpu", "cuda"],
+        help="Device to use for prediction. Default is cpu.",
+    )
+    parser.add_argument(
+        "--work-dir",
+        type=str,
+        default=datetime.now().strftime("%Y-%m-%d-%H-%M-%S") + "-" + str(uuid.uuid4()),
+        help="Working directory for the calculations. "
+        "Defaults to a UUID with timestamp when not set.",
+    )
+    parser.add_argument(
+        "--save-csv",
+        type=str,
+        default="results.csv.gz",
+        help="Save the results to a CSV file. "
+        "Defaults to `results.csv.gz` when not set.",
+    )
+
+
+def add_relax_args(parser: argparse.ArgumentParser):
+    parser.add_argument(
+        "--optimizer",
+        type=str,
+        default="FIRE",
+        help="The optimizer to use. Default is FIRE.",
+    )
+    parser.add_argument(
+        "--filter",
+        type=str,
+        default=None,
+        help="The filter to use.",
+    )
+    parser.add_argument(
+        "--constrain-symmetry",
+        action="store_true",
+        help="Constrain symmetry.",
+    )
+    parser.add_argument(
+        "--fix-axis",
+        type=bool,
+        default=False,
+        nargs="+",
+        help="Fix the axis.",
+    )
+    parser.add_argument(
+        "--pressure-in-GPa",
+        type=float,
+        default=None,
+        help="Pressure in GPa to use for relaxation.",
+    )
+    parser.add_argument(
+        "--fmax",
+        type=float,
+        default=0.01,
+        help="Maximum force tolerance for relaxation.",
+    )
+    parser.add_argument(
+        "--steps",
+        type=int,
+        default=500,
+        help="Maximum number of steps for relaxation.",
+    )
+
+
+def add_phonon_args(parser: argparse.ArgumentParser):
+    parser.add_argument(
+        "--find-prim",
+        action="store_true",
+        help="If find the primitive cell and use it to calculate phonon.",
+    )
+    parser.add_argument(
+        "--amplitude",
+        type=float,
+        default=0.01,
+        help="Magnitude of the finite difference to displace in "
+        "force constant calculation, in Angstrom.",
+    )
+    parser.add_argument(
+        "--supercell-matrix",
+        type=int,
+        nargs=3,
+        default=None,
+        help="Supercell matrix for construct supercell, must be a list of 3 integers.",
+    )
+    parser.add_argument(
+        "--qpoints-mesh",
+        type=int,
+        nargs=3,
+        default=None,
+        help="Qpoint mesh for IBZ integral, must be a list of 3 integers.",
+    )
+    parser.add_argument(
+        "--max-atoms",
+        type=int,
+        default=None,
+        help="Maximum atoms number limitation for the supercell generation.",
+    )
+    parser.add_argument(
+        "--enable-relax",
+        action="store_true",
+        help="Whether to relax the structure before predicting phonon properties.",
+    )
+
+
+def add_moldyn_args(parser: argparse.ArgumentParser):
+    parser.add_argument(
+        "--temperature",
+        type=float,
+        default=300,
+        help="Temperature in Kelvin.",
+    )
+    parser.add_argument(
+        "--timestep",
+        type=float,
+        default=1,
+        help="Timestep in femtoseconds.",
+    )
+    parser.add_argument(
+        "--steps",
+        type=int,
+        default=1000,
+        help="Number of steps for the molecular dynamics simulation.",
+    )
+    parser.add_argument(
+        "--ensemble",
+        type=str,
+        choices=["nvt_berendsen", "nvt_nose_hoover"],
+        default="nvt_nose_hoover",
+        help="Simulation ensemble to use.",
+    )
+    parser.add_argument(
+        "--logfile",
+        type=str,
+        default="md.log",
+        help="Logfile to write the output to. Default is stdout.",
+    )
+    parser.add_argument(
+        "--loginterval",
+        type=int,
+        default=10,
+        help="Log interval for writing the output.",
+    )
+    parser.add_argument(
+        "--trajectory",
+        type=str,
+        default="md.traj",
+        help="Path to the trajectory file.",
+    )
+
+
+def parse_atoms_list(
+    structure_file_list: Union[str, List[str]],
+    mattersim_model: str,
+    device: str = "cpu",
+) -> List[Atoms]:
+    if isinstance(structure_file_list, str):
+        structure_file_list = [structure_file_list]
+
+    calc = MatterSimCalculator(load_path=mattersim_model, device=device)
+    atoms_list = []
+    for structure_file in structure_file_list:
+        atoms_list += ase_read(structure_file, index=":")
+    for atoms in atoms_list:
+        atoms.calc = calc
+    return atoms_list
+
+
+def main():
+    argparser = argparse.ArgumentParser(description="CLI for MatterSim.")
+    subparsers = argparser.add_subparsers(
+        title="Subcommands",
+        description="Valid subcommands",
+        help="Available subcommands",
+    )
+
+    # Sub-command for single-point prediction
+    singlepoint_parser = subparsers.add_parser(
+        "singlepoint", help="Predict single point properties for a list of atoms."
+    )
+    add_common_args(singlepoint_parser)
+    singlepoint_parser.set_defaults(func=singlepoint_cli)
+
+    # Sub-command for relax
+    relax_parser = subparsers.add_parser(
+        "relax", help="Relax a list of atoms structures."
+    )
+    add_common_args(relax_parser)
+    add_relax_args(relax_parser)
+    relax_parser.set_defaults(func=relax_cli)
+
+    # Sub-command for phonon
+    phonon_parser = subparsers.add_parser(
+        "phonon",
+        help="Predict phonon properties for a list of structures.",
+    )
+    add_common_args(phonon_parser)
+    add_relax_args(phonon_parser)
+    add_phonon_args(phonon_parser)
+    phonon_parser.set_defaults(func=phonon_cli)
+
+    # Sub-command for molecular dynamics
+    moldyn_parser = subparsers.add_parser(
+        "moldyn",
+        help="Perform molecular dynamics simulation for a list of structures.",
+    )
+    add_common_args(moldyn_parser)
+    add_moldyn_args(moldyn_parser)
+    moldyn_parser.set_defaults(func=moldyn_cli)
+
+    # Parse arguments
+    args = argparser.parse_args()
+    print(args)
+
+    # Call the function associated with the sub-command
+    if hasattr(args, "func"):
+        args.func(args)
+    else:
+        argparser.print_help()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/data/mp-149_Si2.cif b/tests/data/mp-149_Si2.cif
new file mode 100644
index 0000000..12301a9
--- /dev/null
+++ b/tests/data/mp-149_Si2.cif
@@ -0,0 +1,27 @@
+data_image0
+_chemical_formula_structural       Si2
+_chemical_formula_sum              "Si2"
+_cell_length_a       3.8492784033699095
+_cell_length_b       3.8492794116013456
+_cell_length_c       3.849278
+_cell_angle_alpha    60.00001213094421
+_cell_angle_beta     60.00000346645984
+_cell_angle_gamma    60.00001097545789
+
+_space_group_name_H-M_alt    "P 1"
+_space_group_IT_number       1
+
+loop_
+  _space_group_symop_operation_xyz
+  'x, y, z'
+
+loop_
+  _atom_site_type_symbol
+  _atom_site_label
+  _atom_site_symmetry_multiplicity
+  _atom_site_fract_x
+  _atom_site_fract_y
+  _atom_site_fract_z
+  _atom_site_occupancy
+  Si  Si1       1.0  0.875  0.8749999999999999  0.8750000000000001  1.0000
+  Si  Si2       1.0  0.12500000000000003  0.125  0.12499999999999997  1.0000
diff --git a/tests/data/mp-2998_BaTiO3.cif b/tests/data/mp-2998_BaTiO3.cif
new file mode 100644
index 0000000..4b031b2
--- /dev/null
+++ b/tests/data/mp-2998_BaTiO3.cif
@@ -0,0 +1,35 @@
+data_image0
+_chemical_formula_structural       Ba2Ti2O6
+_chemical_formula_sum              "Ba2 Ti2 O6"
+_cell_length_a       5.666408035279349
+_cell_length_b       5.6664076443956235
+_cell_length_c       5.668373884507181
+_cell_angle_alpha    119.98853562968743
+_cell_angle_beta     119.98853684518572
+_cell_angle_gamma    90.00000405170597
+
+_space_group_name_H-M_alt    "P 1"
+_space_group_IT_number       1
+
+loop_
+  _space_group_symop_operation_xyz
+  'x, y, z'
+
+loop_
+  _atom_site_type_symbol
+  _atom_site_label
+  _atom_site_symmetry_multiplicity
+  _atom_site_fract_x
+  _atom_site_fract_y
+  _atom_site_fract_z
+  _atom_site_occupancy
+  Ba  Ba1       1.0  0.7499999999999999  0.24999999999999997  0.5  1.0000
+  Ba  Ba2       1.0  0.24999999999999994  0.7499999999999999  0.49999999999999994  1.0000
+  Ti  Ti1       1.0  0.5  0.49999999999999994  0.0  1.0000
+  Ti  Ti2       1.0  0.0  0.0  0.0  1.0000
+  O   O1        1.0  0.75046047  0.25046046999999994  0.0  1.0000
+  O   O2        1.0  0.7495395300000001  0.7504604699999999  0.0  1.0000
+  O   O3        1.0  0.25046047  0.24953952999999995  0.0  1.0000
+  O   O4        1.0  0.24953952999999995  0.74953953  0.9999999999999999  1.0000
+  O   O5        1.0  0.25  0.24999999999999997  0.5  1.0000
+  O   O6        1.0  0.75  0.7499999999999999  0.5  1.0000
diff --git a/tests/data/mp-66_C2.cif b/tests/data/mp-66_C2.cif
new file mode 100644
index 0000000..88b7f70
--- /dev/null
+++ b/tests/data/mp-66_C2.cif
@@ -0,0 +1,27 @@
+data_image0
+_chemical_formula_structural       C2
+_chemical_formula_sum              "C2"
+_cell_length_a       2.5178271667719256
+_cell_length_b       2.5178269134954783
+_cell_length_c       2.51782692
+_cell_angle_alpha    59.99999991454223
+_cell_angle_beta     60.00000324214051
+_cell_angle_gamma    59.9999988338285
+
+_space_group_name_H-M_alt    "P 1"
+_space_group_IT_number       1
+
+loop_
+  _space_group_symop_operation_xyz
+  'x, y, z'
+
+loop_
+  _atom_site_type_symbol
+  _atom_site_label
+  _atom_site_symmetry_multiplicity
+  _atom_site_fract_x
+  _atom_site_fract_y
+  _atom_site_fract_z
+  _atom_site_occupancy
+  C   C1        1.0  0.8750000000000001  0.875  0.8749999999999998  1.0000
+  C   C2        1.0  0.125  0.125  0.125  1.0000