diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index 1dd682776..120a36adc 100644
--- a/.github/workflows/benchmark.yml
+++ b/.github/workflows/benchmark.yml
@@ -13,7 +13,12 @@ jobs:
       uses: actions/setup-python@v5
       with:
         python-version: 3.12
-    - run: curl -LsSf https://astral.sh/uv/install.sh | sh
+    - uses: astral-sh/setup-uv@v5
+      with:
+        enable-cache: true
+        cache-dependency-glob: |
+          **/requirements*.txt
+          **/pyproject.toml
     - name: Install dependencies
       run: uv pip install --system .[test,amber,ase,pymatgen,benchmark] rdkit openbabel-wheel
     - name: Run benchmarks
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index a6d23270c..caf4a59c6 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -6,10 +6,10 @@ on:
 
 jobs:
   build:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-22.04
     strategy:
       matrix:
-        python-version: ["3.7", "3.8", "3.12"]
+        python-version: ["3.8", "3.12"]
 
     steps:
     - uses: actions/checkout@v4
@@ -18,13 +18,19 @@ jobs:
       uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
-    - run: curl -LsSf https://astral.sh/uv/install.sh | sh
+    - uses: astral-sh/setup-uv@v5
+      with:
+        enable-cache: true
+        cache-dependency-glob: |
+          **/requirements*.txt
+          **/pyproject.toml
+        cache-suffix: "py${{ matrix.python-version }}"
     - name: Install dependencies
       run: uv pip install --system .[test,amber,ase,pymatgen] coverage ./tests/plugin rdkit openbabel-wheel
     - name: Test
       run: cd tests && coverage run --source=../dpdata -m unittest && cd .. && coverage combine tests/.coverage && coverage report
     - name: Run codecov
-      uses: codecov/codecov-action@v4
+      uses: codecov/codecov-action@v5
       env:
         CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
   pass:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 8da028de5..75f26cf19 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -2,7 +2,7 @@
 # See https://pre-commit.com/hooks.html for more hooks
 repos:
 -   repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.6.0
+    rev: v5.0.0
     hooks:
     # there are many log files in tests
     # TODO: seperate py files and log files
@@ -21,7 +21,7 @@ repos:
 # Python
 -   repo: https://github.com/astral-sh/ruff-pre-commit
     # Ruff version.
-    rev: v0.6.5
+    rev: v0.9.2
     hooks:
     - id: ruff
       args: ["--fix"]
@@ -36,7 +36,7 @@ repos:
       args: ["--write"]
 # Python inside docs
 -   repo: https://github.com/asottile/blacken-docs
-    rev: 1.18.0
+    rev: 1.19.1
     hooks:
     -   id: blacken-docs
 ci:
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 83a360327..79c271f43 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -3,7 +3,7 @@ version: 2
 build:
   os: ubuntu-22.04
   tools:
-    python: "mambaforge-22.9"
+    python: "mambaforge-23.11"
 conda:
   environment: docs/rtd_environment.yml
 sphinx:
diff --git a/README.md b/README.md
index 6f9246b1c..3ca2355c6 100644
--- a/README.md
+++ b/README.md
@@ -8,7 +8,7 @@
 
 ## Installation
 
-DP-GEN only supports Python 3.7 and above. You can [setup a conda/pip environment](https://docs.deepmodeling.com/faq/conda.html), and then use one of the following methods to install DP-GEN:
+dpdata only supports Python 3.8 and above. You can [setup a conda/pip environment](https://docs.deepmodeling.com/faq/conda.html), and then use one of the following methods to install dpdata:
 
 - Install via pip: `pip install dpdata`
 - Install via conda: `conda install -c conda-forge dpdata`
diff --git a/docs/conf.py b/docs/conf.py
index a21d19d56..263cb5507 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -24,7 +24,7 @@
 # -- Project information -----------------------------------------------------
 
 project = "dpdata"
-copyright = "2019-%d, DeepModeling " % date.today().year
+copyright = "2019-%d, DeepModeling " % date.today().year  # noqa: UP031
 author = "Han Wang"
 
 # The short X.Y version
@@ -44,7 +44,7 @@
 # ones.
 extensions = [
     "deepmodeling_sphinx",
-    "sphinx_rtd_theme",
+    "sphinx_book_theme",
     "sphinx.ext.mathjax",
     "sphinx.ext.viewcode",
     "sphinx.ext.intersphinx",
@@ -87,7 +87,7 @@
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = "sphinx_rtd_theme"
+html_theme = "sphinx_book_theme"
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
diff --git a/docs/installation.md b/docs/installation.md
index 064f91331..1b24e0051 100644
--- a/docs/installation.md
+++ b/docs/installation.md
@@ -1,6 +1,6 @@
 # Installation
 
-DP-GEN only supports Python 3.7 and above. You can [setup a conda/pip environment](https://docs.deepmodeling.com/faq/conda.html), and then use one of the following methods to install DP-GEN:
+dpdata only supports Python 3.8 and above. You can [setup a conda/pip environment](https://docs.deepmodeling.com/faq/conda.html), and then use one of the following methods to install dpdata:
 
 - Install via pip: `pip install dpdata`
 - Install via conda: `conda install -c conda-forge dpdata`
diff --git a/docs/make_format.py b/docs/make_format.py
index 2b3c03c67..e9c1f60d3 100644
--- a/docs/make_format.py
+++ b/docs/make_format.py
@@ -2,14 +2,9 @@
 
 import csv
 import os
-import sys
 from collections import defaultdict
 from inspect import Parameter, Signature, cleandoc, signature
-
-if sys.version_info >= (3, 8):
-    from typing import Literal
-else:
-    from typing_extensions import Literal
+from typing import Literal
 
 from numpydoc.docscrape import Parameter as numpydoc_Parameter
 from numpydoc.docscrape_sphinx import SphinxDocString
diff --git a/docs/rtd_environment.yml b/docs/rtd_environment.yml
index 7eae4aad6..50572dca7 100644
--- a/docs/rtd_environment.yml
+++ b/docs/rtd_environment.yml
@@ -2,6 +2,7 @@ name: dpdata
 channels:
   - conda-forge
 dependencies:
-  - mamba
+  - python <3.13
+  - mamba <2
   - pip:
     - ..[docs]
diff --git a/dpdata/abacus/md.py b/dpdata/abacus/md.py
index bbb023439..dbe86735f 100644
--- a/dpdata/abacus/md.py
+++ b/dpdata/abacus/md.py
@@ -56,9 +56,9 @@ def get_coords_from_dump(dumplines, natoms):
     if "VIRIAL" in dumplines[6]:
         calc_stress = True
         check_line = 10
-    assert (
-        "POSITION" in dumplines[check_line]
-    ), "keywords 'POSITION' cannot be found in the 6th line. Please check."
+    assert "POSITION" in dumplines[check_line], (
+        "keywords 'POSITION' cannot be found in the 6th line. Please check."
+    )
     if "FORCE" in dumplines[check_line]:
         calc_force = True
 
@@ -68,7 +68,7 @@ def get_coords_from_dump(dumplines, natoms):
     else:
         nframes_dump = int(nlines / (total_natoms + 9))
     assert nframes_dump > 0, (
-        "Number of lines in MD_dump file = %d. Number of atoms = %d. The MD_dump file is incomplete."
+        "Number of lines in MD_dump file = %d. Number of atoms = %d. The MD_dump file is incomplete."  # noqa: UP031
         % (nlines, total_natoms)
     )
     cells = np.zeros([nframes_dump, 3, 3])
@@ -125,7 +125,7 @@ def get_coords_from_dump(dumplines, natoms):
                     )
             iframe += 1
     assert iframe == nframes_dump, (
-        "iframe=%d, nframe_dump=%d. Number of frames does not match number of lines in MD_dump."
+        "iframe=%d, nframe_dump=%d. Number of frames does not match number of lines in MD_dump."  # noqa: UP031
         % (iframe, nframes_dump)
     )
     stresses *= kbar2evperang3
@@ -145,7 +145,7 @@ def get_energy(outlines, ndump, dump_freq):
                 energy.append(np.nan)
             nenergy += 1
     assert ndump == len(energy), (
-        "Number of total energies in running_md.log = %d. Number of frames in MD_dump = %d. Please check."
+        "Number of total energies in running_md.log = %d. Number of frames in MD_dump = %d. Please check."  # noqa: UP031
         % (len(energy), ndump)
     )
     energy = np.array(energy)
@@ -167,7 +167,7 @@ def get_frame(fname):
     with open_file(geometry_path_in) as fp:
         geometry_inlines = fp.read().split("\n")
     celldm, cell = get_cell(geometry_inlines)
-    atom_names, natoms, types, coords = get_coords(
+    atom_names, natoms, types, coords, move, magmom = get_coords(
         celldm, cell, geometry_inlines, inlines
     )
     # This coords is not to be used.
@@ -191,7 +191,7 @@ def get_frame(fname):
             force = np.delete(force, i - ndump, axis=0)
             stress = np.delete(stress, i - ndump, axis=0)
             energy = np.delete(energy, i - ndump, axis=0)
-            unconv_stru += "%d " % i
+            unconv_stru += "%d " % i  # noqa: UP031
     ndump = len(energy)
     if unconv_stru != "":
         warnings.warn(f"Structure {unconv_stru} are unconverged and not collected!")
@@ -220,6 +220,8 @@ def get_frame(fname):
     if len(magmom) > 0:
         data["spins"] = magmom
     if len(magforce) > 0:
-        data["mag_forces"] = magforce
+        data["force_mags"] = magforce
+    if len(move) > 0:
+        data["move"] = move
 
     return data
diff --git a/dpdata/abacus/relax.py b/dpdata/abacus/relax.py
index be80bc90a..8f90040ef 100644
--- a/dpdata/abacus/relax.py
+++ b/dpdata/abacus/relax.py
@@ -47,7 +47,7 @@ def get_coords_from_log(loglines, natoms):
             natoms_log += int(line.split()[-1])
 
     assert natoms_log > 0 and natoms_log == natoms, (
-        "ERROR: detected atom number in log file is %d" % natoms
+        "ERROR: detected atom number in log file is %d" % natoms  # noqa: UP031
     )
 
     energy = []
@@ -76,7 +76,7 @@ def get_coords_from_log(loglines, natoms):
                         list(map(lambda x: float(x) * a0, loglines[i + k].split()[1:4]))
                     )
             else:
-                assert False, "Unrecongnized coordinate type, %s, line:%d" % (
+                assert False, "Unrecongnized coordinate type, %s, line:%d" % (  # noqa: UP031
                     loglines[i].split()[0],
                     i,
                 )
@@ -183,7 +183,7 @@ def get_frame(fname):
     with open_file(geometry_path_in) as fp:
         geometry_inlines = fp.read().split("\n")
     celldm, cell = get_cell(geometry_inlines)
-    atom_names, natoms, types, coord_tmp = get_coords(
+    atom_names, natoms, types, coord_tmp, move, magmom = get_coords(
         celldm, cell, geometry_inlines, inlines
     )
 
@@ -217,6 +217,8 @@ def get_frame(fname):
     if len(magmom) > 0:
         data["spins"] = magmom
     if len(magforce) > 0:
-        data["mag_forces"] = magforce
+        data["force_mags"] = magforce
+    if len(move) > 0:
+        data["move"] = move
 
     return data
diff --git a/dpdata/abacus/scf.py b/dpdata/abacus/scf.py
index 9919e9128..5c1f283cc 100644
--- a/dpdata/abacus/scf.py
+++ b/dpdata/abacus/scf.py
@@ -258,6 +258,51 @@ def parse_stru_pos(pos_line):
     return pos, move, velocity, magmom, angle1, angle2, constrain, lambda1
 
 
+def get_atom_mag_cartesian(atommag, angle1, angle2):
+    """Transform atommag, angle1, angle2 to magmom in cartesian coordinates.
+
+    Parameters
+    ----------
+    atommag : float/list of float/None
+        Atom magnetic moment.
+    angle1 : float/None
+        value of angle1.
+    angle2 : float/None
+        value of angle2.
+    ABACUS support defining mag, angle1, angle2 at the same time.
+    angle1 is the angle between z-axis and real spin (in degrees).
+    angle2 is the angle between x-axis and real spin projection in xy-plane (in degrees).
+    If only mag is defined, then transfer it to magmom directly.
+    And if mag, angle1, angle2 are defined, then mag is only the norm of magmom, and the direction is defined by angle1 and angle2.
+    """
+    if atommag is None:
+        return None
+    if not (isinstance(atommag, list) or isinstance(atommag, float)):
+        raise RuntimeError(f"Invalid atommag: {atommag}")
+
+    if angle1 is None and angle2 is None:
+        if isinstance(atommag, list):
+            return atommag
+        else:
+            return [0, 0, atommag]
+    else:
+        a1 = 0
+        a2 = 0
+        if angle1 is not None:
+            a1 = angle1
+        if angle2 is not None:
+            a2 = angle2
+        if isinstance(atommag, list):
+            mag_norm = np.linalg.norm(atommag)
+        else:
+            mag_norm = atommag
+        return [
+            mag_norm * np.sin(np.radians(a1)) * np.cos(np.radians(a2)),
+            mag_norm * np.sin(np.radians(a1)) * np.sin(np.radians(a2)),
+            mag_norm * np.cos(np.radians(a1)),
+        ]
+
+
 def get_coords(celldm, cell, geometry_inlines, inlines=None):
     coords_lines = get_stru_block(geometry_inlines, "ATOMIC_POSITIONS")
     # assuming that ATOMIC_POSITIONS is at the bottom of the STRU file
@@ -268,9 +313,7 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None):
     coords = []  # coordinations of atoms
     move = []  # move flag of each atom
     velocity = []  # velocity of each atom
-    mag = []  # magnetic moment of each atom
-    angle1 = []  # angle1 of each atom
-    angle2 = []  # angle2 of each atom
+    mags = []  # magnetic moment of each atom
     sc = []  # spin constraint flag of each atom
     lambda_ = []  # lambda of each atom
 
@@ -278,6 +321,7 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None):
     line_idx = 1  # starting line of first element
     for it in range(ntype):
         atom_names.append(coords_lines[line_idx].split()[0])
+        atom_type_mag = float(coords_lines[line_idx + 1].split()[0])
         line_idx += 2
         atom_numbs.append(int(coords_lines[line_idx].split()[0]))
         line_idx += 1
@@ -299,18 +343,23 @@ def get_coords(celldm, cell, geometry_inlines, inlines=None):
             coords.append(xyz)
             atom_types.append(it)
 
-            move.append(imove)
+            if imove is not None:
+                move.append(imove)
             velocity.append(ivelocity)
-            mag.append(imagmom)
-            angle1.append(iangle1)
-            angle2.append(iangle2)
             sc.append(iconstrain)
             lambda_.append(ilambda1)
 
+            # calculate the magnetic moment in cartesian coordinates
+            mag = get_atom_mag_cartesian(imagmom, iangle1, iangle2)
+            if mag is None:
+                mag = [0, 0, atom_type_mag]
+            mags.append(mag)
+
             line_idx += 1
     coords = np.array(coords)  # need transformation!!!
     atom_types = np.array(atom_types)
-    return atom_names, atom_numbs, atom_types, coords
+    move = np.array(move, dtype=bool)
+    return atom_names, atom_numbs, atom_types, coords, move, mags
 
 
 def get_energy(outlines):
@@ -429,14 +478,20 @@ def get_mag_force(outlines):
             j = i + 2
             mag = []
             while "-------------------------" not in outlines[j]:
-                mag.append([float(ii) for ii in outlines[j].split()[1:]])
+                imag = [float(ii) for ii in outlines[j].split()[1:]]
+                if len(imag) == 1:
+                    imag = [0, 0, imag[0]]
+                mag.append(imag)
                 j += 1
             mags.append(mag)
         if "Magnetic force (eV/uB)" in line:
             j = i + 2
             magforce = []
             while "-------------------------" not in outlines[j]:
-                magforce.append([float(ii) for ii in outlines[j].split()[1:]])
+                imagforce = [float(ii) for ii in outlines[j].split()[1:]]
+                if len(imagforce) == 1:
+                    imagforce = [0, 0, imagforce[0]]
+                magforce.append(imagforce)
                 j += 1
             magforces.append(magforce)
     return np.array(mags), np.array(magforces)
@@ -477,9 +532,12 @@ def get_frame(fname):
         outlines = fp.read().split("\n")
 
     celldm, cell = get_cell(geometry_inlines)
-    atom_names, natoms, types, coords = get_coords(
-        celldm, cell, geometry_inlines, inlines
+    atom_names, natoms, types, coords, move, magmom = (
+        get_coords(  # here the magmom is the initial magnetic moment in STRU
+            celldm, cell, geometry_inlines, inlines
+        )
     )
+
     magmom, magforce = get_mag_force(outlines)
     if len(magmom) > 0:
         magmom = magmom[-1:]
@@ -509,7 +567,9 @@ def get_frame(fname):
     if len(magmom) > 0:
         data["spins"] = magmom
     if len(magforce) > 0:
-        data["mag_forces"] = magforce
+        data["force_mags"] = magforce
+    if len(move) > 0:
+        data["move"] = move[np.newaxis, :, :]
     # print("atom_names = ", data['atom_names'])
     # print("natoms = ", data['atom_numbs'])
     # print("types = ", data['atom_types'])
@@ -561,7 +621,7 @@ def get_frame_from_stru(fname):
     nele = get_nele_from_stru(geometry_inlines)
     inlines = [f"ntype {nele}"]
     celldm, cell = get_cell(geometry_inlines)
-    atom_names, natoms, types, coords = get_coords(
+    atom_names, natoms, types, coords, move, magmom = get_coords(
         celldm, cell, geometry_inlines, inlines
     )
     data = {}
@@ -571,6 +631,9 @@ def get_frame_from_stru(fname):
     data["cells"] = cell[np.newaxis, :, :]
     data["coords"] = coords[np.newaxis, :, :]
     data["orig"] = np.zeros(3)
+    data["spins"] = np.array([magmom])
+    if len(move) > 0:
+        data["move"] = move[np.newaxis, :, :]
 
     return data
 
@@ -601,7 +664,7 @@ def make_unlabeled_stru(
         System data
     frame_idx : int
         The index of the frame to dump
-    pp_file : list of string or dict, optional
+    pp_file : list of string or dict
         List of pseudo potential files, or a dictionary of pseudo potential files for each atomnames
     numerical_orbital : list of string or dict, optional
         List of orbital files, or a dictionary of orbital files for each atomnames
@@ -609,8 +672,8 @@ def make_unlabeled_stru(
         numerical descriptor file
     mass : list of float, optional
         List of atomic masses
-    move : list of list of bool, optional
-        List of the move flag of each xyz direction of each atom
+    move : list of (list of list of bool), optional
+        List of the move flag of each xyz direction of each atom for each frame
     velocity : list of list of float, optional
         List of the velocity of each xyz direction of each atom
     mag : list of (list of float or float), optional
@@ -628,6 +691,8 @@ def make_unlabeled_stru(
     link_file : bool, optional
         Whether to link the pseudo potential files and orbital files in the STRU file.
         If True, then only filename will be written in the STRU file, and make a soft link to the real file.
+    dest_dir : str, optional
+        The destination directory to make the soft link of the pseudo potential files and orbital files.
     For velocity, mag, angle1, angle2, sc, and lambda_, if the value is None, then the corresponding information will not be written.
     ABACUS support defining "mag" and "angle1"/"angle2" at the same time, and in this case, the "mag" only define the norm of the magnetic moment, and "angle1" and "angle2" define the direction of the magnetic moment.
     If data has spins, then it will be written as mag to STRU file; while if mag is passed at the same time, then mag will be used.
@@ -655,6 +720,23 @@ def ndarray2list(i):
         else:
             return i
 
+    def process_file_input(file_input, atom_names, input_name):
+        # For pp_file and numerical_orbital, process the file input, and return a list of file names
+        # file_input can be a list of file names, or a dictionary of file names for each atom names
+        if isinstance(file_input, (list, tuple)):
+            if len(file_input) != len(atom_names):
+                raise ValueError(
+                    f"{input_name} length is not equal to the number of atom types"
+                )
+            return file_input
+        elif isinstance(file_input, dict):
+            for element in atom_names:
+                if element not in file_input:
+                    raise KeyError(f"{input_name} does not contain {element}")
+            return [file_input[element] for element in atom_names]
+        else:
+            raise ValueError(f"Invalid {input_name}: {file_input}")
+
     if link_file and dest_dir is None:
         print(
             "WARNING: make_unlabeled_stru: link_file is True, but dest_dir is None. Will write the filename to STRU but not making soft link."
@@ -665,6 +747,9 @@ def ndarray2list(i):
     if mag is None and data.get("spins") is not None and len(data["spins"]) > 0:
         mag = data["spins"][frame_idx]
 
+    if move is None and data.get("move", None) is not None and len(data["move"]) > 0:
+        move = data["move"][frame_idx]
+
     atom_numbs = sum(data["atom_numbs"])
     for key in [move, velocity, mag, angle1, angle2, sc, lambda_]:
         if key is not None:
@@ -681,7 +766,15 @@ def ndarray2list(i):
     # ATOMIC_SPECIES block
     out = "ATOMIC_SPECIES\n"
     if pp_file is not None:
-        pp_file = ndarray2list(pp_file)
+        ppfiles = process_file_input(
+            ndarray2list(pp_file), data["atom_names"], "pp_file"
+        )
+    else:
+        warnings.warn(
+            "pp_file is not provided, will use empty string for pseudo potential file."
+        )
+        ppfiles = [""] * len(data["atom_names"])
+
     for iele in range(len(data["atom_names"])):
         if data["atom_numbs"][iele] == 0:
             continue
@@ -690,57 +783,37 @@ def ndarray2list(i):
             out += f"{mass[iele]:.3f} "
         else:
             out += "1 "
-        if pp_file is not None:
-            if isinstance(pp_file, (list, tuple)):
-                ipp_file = pp_file[iele]
-            elif isinstance(pp_file, dict):
-                if data["atom_names"][iele] not in pp_file:
-                    print(
-                        f"ERROR: make_unlabeled_stru: pp_file does not contain {data['atom_names'][iele]}"
-                    )
-                    ipp_file = None
-                else:
-                    ipp_file = pp_file[data["atom_names"][iele]]
-            else:
-                ipp_file = None
-            if ipp_file is not None:
-                if not link_file:
-                    out += ipp_file
-                else:
-                    out += os.path.basename(ipp_file.rstrip("/"))
-                    if dest_dir is not None:
-                        _link_file(dest_dir, ipp_file)
 
+        ipp_file = ppfiles[iele]
+        if ipp_file != "":
+            if not link_file:
+                out += ipp_file
+            else:
+                out += os.path.basename(ipp_file.rstrip("/"))
+                if dest_dir is not None:
+                    _link_file(dest_dir, ipp_file)
         out += "\n"
     out += "\n"
 
     # NUMERICAL_ORBITAL block
     if numerical_orbital is not None:
-        assert len(numerical_orbital) == len(data["atom_names"])
         numerical_orbital = ndarray2list(numerical_orbital)
+        orbfiles = process_file_input(
+            numerical_orbital, data["atom_names"], "numerical_orbital"
+        )
+        orbfiles = [
+            orbfiles[i]
+            for i in range(len(data["atom_names"]))
+            if data["atom_numbs"][i] != 0
+        ]
         out += "NUMERICAL_ORBITAL\n"
-        for iele in range(len(data["atom_names"])):
-            if data["atom_numbs"][iele] == 0:
-                continue
-            if isinstance(numerical_orbital, (list, tuple)):
-                inum_orbital = numerical_orbital[iele]
-            elif isinstance(numerical_orbital, dict):
-                if data["atom_names"][iele] not in numerical_orbital:
-                    print(
-                        f"ERROR: make_unlabeled_stru: numerical_orbital does not contain {data['atom_names'][iele]}"
-                    )
-                    inum_orbital = None
-                else:
-                    inum_orbital = numerical_orbital[data["atom_names"][iele]]
+        for iorb in orbfiles:
+            if not link_file:
+                out += iorb
             else:
-                inum_orbital = None
-            if inum_orbital is not None:
-                if not link_file:
-                    out += inum_orbital
-                else:
-                    out += os.path.basename(inum_orbital.rstrip("/"))
-                    if dest_dir is not None:
-                        _link_file(dest_dir, inum_orbital)
+                out += os.path.basename(iorb.rstrip("/"))
+                if dest_dir is not None:
+                    _link_file(dest_dir, iorb)
             out += "\n"
         out += "\n"
 
diff --git a/dpdata/ase_calculator.py b/dpdata/ase_calculator.py
index 1de760a5a..94f5073cb 100644
--- a/dpdata/ase_calculator.py
+++ b/dpdata/ase_calculator.py
@@ -62,7 +62,8 @@ def calculate(
         self.results["energy"] = data["energies"][0]
         # see https://gitlab.com/ase/ase/-/merge_requests/2485
         self.results["free_energy"] = data["energies"][0]
-        self.results["forces"] = data["forces"][0]
+        if "forces" in data:
+            self.results["forces"] = data["forces"][0]
         if "virials" in data:
             self.results["virial"] = data["virials"][0].reshape(3, 3)
 
diff --git a/dpdata/data_type.py b/dpdata/data_type.py
index 851a65c95..8cea28c18 100644
--- a/dpdata/data_type.py
+++ b/dpdata/data_type.py
@@ -122,7 +122,7 @@ def check(self, system: System):
                 elif isinstance(data, list):
                     if len(shape) and shape[0] != len(data):
                         raise DataError(
-                            "Length of %s is %d, but expected %d"
+                            "Length of %s is %d, but expected %d"  # noqa: UP031
                             % (self.name, len(data), shape[0])
                         )
                 else:
diff --git a/dpdata/deepmd/comp.py b/dpdata/deepmd/comp.py
index 67eddd9f1..410d789e1 100644
--- a/dpdata/deepmd/comp.py
+++ b/dpdata/deepmd/comp.py
@@ -119,7 +119,7 @@ def dump(folder, data, set_size=5000, comp_prec=np.float32, remove_sets=True):
     for ii in range(nsets):
         set_stt = ii * set_size
         set_end = (ii + 1) * set_size
-        set_folder = os.path.join(folder, "set.%03d" % ii)
+        set_folder = os.path.join(folder, "set.%03d" % ii)  # noqa: UP031
         os.makedirs(set_folder)
     try:
         os.remove(os.path.join(folder, "nopbc"))
@@ -158,5 +158,5 @@ def dump(folder, data, set_size=5000, comp_prec=np.float32, remove_sets=True):
         for ii in range(nsets):
             set_stt = ii * set_size
             set_end = (ii + 1) * set_size
-            set_folder = os.path.join(folder, "set.%03d" % ii)
+            set_folder = os.path.join(folder, "set.%03d" % ii)  # noqa: UP031
             np.save(os.path.join(set_folder, dtype.deepmd_name), ddata[set_stt:set_end])
diff --git a/dpdata/deepmd/hdf5.py b/dpdata/deepmd/hdf5.py
index cbdaa632b..c2b3bd424 100644
--- a/dpdata/deepmd/hdf5.py
+++ b/dpdata/deepmd/hdf5.py
@@ -59,7 +59,7 @@ def to_system_data(
     else:
         my_type_map = []
         for ii in range(ntypes):
-            my_type_map.append("Type_%d" % ii)
+            my_type_map.append("Type_%d" % ii)  # noqa: UP031
     assert len(my_type_map) >= len(data["atom_numbs"])
     for ii in range(len(data["atom_numbs"])):
         data["atom_names"].append(my_type_map[ii])
@@ -217,7 +217,7 @@ def dump(
     for ii in range(nsets):
         set_stt = ii * set_size
         set_end = (ii + 1) * set_size
-        set_folder = g.create_group("set.%03d" % ii)
+        set_folder = g.create_group("set.%03d" % ii)  # noqa: UP031
         for dt, prop in data_types.items():
             if dt in reshaped_data:
                 set_folder.create_dataset(
diff --git a/dpdata/deepmd/mixed.py b/dpdata/deepmd/mixed.py
index 2aff3a28a..fe5bdaa86 100644
--- a/dpdata/deepmd/mixed.py
+++ b/dpdata/deepmd/mixed.py
@@ -15,9 +15,9 @@ def to_system_data(folder, type_map=None, labels=True):
     if type_map is not None:
         assert isinstance(type_map, list)
         missing_type = [i for i in old_type_map if i not in type_map]
-        assert (
-            not missing_type
-        ), f"These types are missing in selected type_map: {missing_type} !"
+        assert not missing_type, (
+            f"These types are missing in selected type_map: {missing_type} !"
+        )
         index_map = np.array([type_map.index(i) for i in old_type_map])
         data["atom_names"] = type_map.copy()
     else:
diff --git a/dpdata/deepmd/raw.py b/dpdata/deepmd/raw.py
index 81b01cc06..50dc5afd3 100644
--- a/dpdata/deepmd/raw.py
+++ b/dpdata/deepmd/raw.py
@@ -27,7 +27,7 @@ def load_type(folder, type_map=None):
     else:
         my_type_map = []
         for ii in range(ntypes):
-            my_type_map.append("Type_%d" % ii)
+            my_type_map.append("Type_%d" % ii)  # noqa: UP031
     data["atom_names"] = my_type_map
     data["atom_numbs"] = []
     for ii, _ in enumerate(data["atom_names"]):
diff --git a/dpdata/driver.py b/dpdata/driver.py
index b5ff53403..b63c417af 100644
--- a/dpdata/driver.py
+++ b/dpdata/driver.py
@@ -166,7 +166,8 @@ def label(self, data: dict) -> dict:
                 labeled_data = lb_data.copy()
             else:
                 labeled_data["energies"] += lb_data["energies"]
-                labeled_data["forces"] += lb_data["forces"]
+                if "forces" in labeled_data and "forces" in lb_data:
+                    labeled_data["forces"] += lb_data["forces"]
                 if "virials" in labeled_data and "virials" in lb_data:
                     labeled_data["virials"] += lb_data["virials"]
         return labeled_data
diff --git a/dpdata/gaussian/gjf.py b/dpdata/gaussian/gjf.py
index b83dad1c2..419ec354c 100644
--- a/dpdata/gaussian/gjf.py
+++ b/dpdata/gaussian/gjf.py
@@ -186,8 +186,8 @@ def make_gaussian_input(
             mult_frags.append(detect_multiplicity(np.array(symbols)[idx]))
         if use_fragment_guesses:
             multiplicity = sum(mult_frags) - frag_numb + 1 - charge % 2
-            chargekeywords_frag = "%d %d" % (charge, multiplicity) + "".join(
-                [" %d %d" % (charge, mult_frag) for mult_frag in mult_frags]
+            chargekeywords_frag = "%d %d" % (charge, multiplicity) + "".join(  # noqa: UP031
+                [" %d %d" % (charge, mult_frag) for mult_frag in mult_frags]  # noqa: UP031
             )
         else:
             multi_frags = np.array(mult_frags)
@@ -239,10 +239,10 @@ def make_gaussian_input(
     for ii, (symbol, coordinate) in enumerate(zip(symbols, coordinates)):
         if use_fragment_guesses:
             buff.append(
-                "%s(Fragment=%d) %f %f %f" % (symbol, frag_index[ii] + 1, *coordinate)
+                "%s(Fragment=%d) %f %f %f" % (symbol, frag_index[ii] + 1, *coordinate)  # noqa: UP031
             )
         else:
-            buff.append("{} {:f} {:f} {:f}".format(symbol, *coordinate))
+            buff.append("{} {:f} {:f} {:f}".format(symbol, *coordinate))  # noqa: UP031
     if not sys_data.get("nopbc", False):
         # PBC condition
         cell = sys_data["cells"][0]
diff --git a/dpdata/lammps/dump.py b/dpdata/lammps/dump.py
index 72fee27df..fe549b956 100644
--- a/dpdata/lammps/dump.py
+++ b/dpdata/lammps/dump.py
@@ -197,7 +197,91 @@ def load_file(fname: FileType, begin=0, step=1):
                 buff.append(line)
 
 
-def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False):
+def get_spin_keys(inputfile):
+    """
+    Read input file and get the keys for spin info in dump.
+
+    Parameters
+    ----------
+    inputfile : str
+        Path to the input file.
+
+    Returns
+    -------
+    list or None
+        List of spin info keys if found, None otherwise.
+    """
+    if inputfile is None:
+        return None
+
+    if not os.path.isfile(inputfile):
+        warnings.warn(f"Input file {inputfile} not found.")
+        return None
+
+    with open(inputfile) as f:
+        for line in f.readlines():
+            ls = line.split()
+            if (
+                len(ls) > 7
+                and ls[0] == "compute"
+                and all(key in ls for key in ["sp", "spx", "spy", "spz"])
+            ):
+                compute_name = ls[1]
+                return [
+                    f"c_{compute_name}[{ls.index(key) - 3}]"
+                    for key in ["sp", "spx", "spy", "spz"]
+                ]
+
+    return None
+
+
+def get_spin(lines, spin_keys):
+    """
+    Get the spin info from the dump file.
+
+    Parameters
+    ----------
+    lines : list
+        The content of the dump file.
+    spin_keys : list
+        The keys for spin info in dump file.
+    the spin info is stored in sp, spx, spy, spz or spin_keys, which is the spin norm and the spin vector
+    1 1 0.00141160 5.64868599 0.01005602 1.54706291 0.00000000 0.00000000 1.00000000 -1.40772100 -2.03739417 -1522.64797384 -0.00397809 -0.00190426 -0.00743976
+    """
+    blk, head = _get_block(lines, "ATOMS")
+    heads = head.split()
+
+    if spin_keys is not None and all(i in heads for i in spin_keys):
+        key = spin_keys
+    else:
+        return None
+
+    try:
+        idx_id = heads.index("id") - 2
+        idx_sp, idx_spx, idx_spy, idx_spz = (heads.index(k) - 2 for k in key)
+
+        norm = []
+        vec = []
+        atom_ids = []
+        for line in blk:
+            words = line.split()
+            norm.append([float(words[idx_sp])])
+            vec.append(
+                [float(words[idx_spx]), float(words[idx_spy]), float(words[idx_spz])]
+            )
+            atom_ids.append(int(words[idx_id]))
+
+        spin = np.array(norm) * np.array(vec)
+        atom_ids, spin = zip(*sorted(zip(atom_ids, spin)))
+        return np.array(spin)
+    except (ValueError, IndexError) as e:
+        warnings.warn(f"Error processing spin data: {str(e)}")
+        return None
+
+
+def system_data(
+    lines, type_map=None, type_idx_zero=True, unwrap=False, input_file=None
+):
     array_lines = split_traj(lines)
     lines = array_lines[0]
     system = {}
@@ -205,7 +289,7 @@ def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False):
     system["atom_names"] = []
     if type_map is None:
         for ii in range(len(system["atom_numbs"])):
-            system["atom_names"].append("TYPE_%d" % ii)
+            system["atom_names"].append("TYPE_%d" % ii)  # noqa: UP031
     else:
         assert len(type_map) >= len(system["atom_numbs"])
         for ii in range(len(system["atom_numbs"])):
@@ -216,6 +300,12 @@ def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False):
     system["cells"] = [np.array(cell)]
     system["atom_types"] = get_atype(lines, type_idx_zero=type_idx_zero)
     system["coords"] = [safe_get_posi(lines, cell, np.array(orig), unwrap)]
+    spin_keys = get_spin_keys(input_file)
+    spin = get_spin(lines, spin_keys)
+    has_spin = False
+    if spin is not None:
+        system["spins"] = [spin]
+        has_spin = True
     for ii in range(1, len(array_lines)):
         bounds, tilt = get_dumpbox(array_lines[ii])
         orig, cell = dumpbox2box(bounds, tilt)
@@ -228,6 +318,18 @@ def system_data(lines, type_map=None, type_idx_zero=True, unwrap=False):
         system["coords"].append(
             safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap)[idx]
         )
+        if has_spin:
+            spin = get_spin(array_lines[ii], spin_keys)
+            if spin is not None:
+                system["spins"].append(spin[idx])
+            else:
+                warnings.warn(
+                    f"Warning: spin info is not found in frame {ii}, remove spin info."
+                )
+                system.pop("spins")
+                has_spin = False
+    if has_spin:
+        system["spins"] = np.array(system["spins"])
     system["cells"] = np.array(system["cells"])
     system["coords"] = np.array(system["coords"])
     return system
diff --git a/dpdata/lammps/lmp.py b/dpdata/lammps/lmp.py
index 604b18d12..d9f2f1d51 100644
--- a/dpdata/lammps/lmp.py
+++ b/dpdata/lammps/lmp.py
@@ -127,6 +127,19 @@ def get_posi(lines):
     return np.array(posis)
 
 
+def get_spins(lines):
+    atom_lines = get_atoms(lines)
+    if len(atom_lines[0].split()) < 8:
+        return None
+    spins_ori = []
+    spins_norm = []
+    for ii in atom_lines:
+        iis = ii.split()
+        spins_ori.append([float(jj) for jj in iis[5:8]])
+        spins_norm.append([float(iis[-1])])
+    return np.array(spins_ori) * np.array(spins_norm)
+
+
 def get_lmpbox(lines):
     box_info = []
     tilt = np.zeros(3)
@@ -154,7 +167,7 @@ def system_data(lines, type_map=None, type_idx_zero=True):
     system["atom_names"] = []
     if type_map is None:
         for ii in range(len(system["atom_numbs"])):
-            system["atom_names"].append("Type_%d" % ii)
+            system["atom_names"].append("Type_%d" % ii)  # noqa: UP031
     else:
         assert len(type_map) >= len(system["atom_numbs"])
         for ii in range(len(system["atom_numbs"])):
@@ -168,6 +181,11 @@ def system_data(lines, type_map=None, type_idx_zero=True):
     system["coords"] = [get_posi(lines)]
     system["cells"] = np.array(system["cells"])
     system["coords"] = np.array(system["coords"])
+
+    spins = get_spins(lines)
+    if spins is not None:
+        system["spins"] = np.array([spins])
+
     return system
 
 
@@ -180,27 +198,27 @@ def from_system_data(system, f_idx=0):
     ret += "\n"
     natoms = sum(system["atom_numbs"])
     ntypes = len(system["atom_numbs"])
-    ret += "%d atoms\n" % natoms
-    ret += "%d atom types\n" % ntypes
+    ret += "%d atoms\n" % natoms  # noqa: UP031
+    ret += "%d atom types\n" % ntypes  # noqa: UP031
     ret += (ptr_float_fmt + " " + ptr_float_fmt + " xlo xhi\n") % (
         0,
         system["cells"][f_idx][0][0],
-    )
+    )  # noqa: UP031
     ret += (ptr_float_fmt + " " + ptr_float_fmt + " ylo yhi\n") % (
         0,
         system["cells"][f_idx][1][1],
-    )
+    )  # noqa: UP031
     ret += (ptr_float_fmt + " " + ptr_float_fmt + " zlo zhi\n") % (
         0,
         system["cells"][f_idx][2][2],
-    )
+    )  # noqa: UP031
     ret += (
         ptr_float_fmt + " " + ptr_float_fmt + " " + ptr_float_fmt + " xy xz yz\n"
     ) % (
         system["cells"][f_idx][1][0],
         system["cells"][f_idx][2][0],
         system["cells"][f_idx][2][1],
-    )
+    )  # noqa: UP031
     ret += "\n"
     ret += "Atoms # atomic\n"
     ret += "\n"
@@ -215,15 +233,56 @@ def from_system_data(system, f_idx=0):
         + " "
         + ptr_float_fmt
         + "\n"
-    )
+    )  # noqa: UP031
+
+    if "spins" in system:
+        coord_fmt = (
+            coord_fmt.strip("\n")
+            + " "
+            + ptr_float_fmt
+            + " "
+            + ptr_float_fmt
+            + " "
+            + ptr_float_fmt
+            + " "
+            + ptr_float_fmt
+            + "\n"
+        )  # noqa: UP031
+        spins_norm = np.linalg.norm(system["spins"][f_idx], axis=1)
     for ii in range(natoms):
-        ret += coord_fmt % (
-            ii + 1,
-            system["atom_types"][ii] + 1,
-            system["coords"][f_idx][ii][0] - system["orig"][0],
-            system["coords"][f_idx][ii][1] - system["orig"][1],
-            system["coords"][f_idx][ii][2] - system["orig"][2],
-        )
+        if "spins" in system:
+            if spins_norm[ii] != 0:
+                ret += coord_fmt % (
+                    ii + 1,
+                    system["atom_types"][ii] + 1,
+                    system["coords"][f_idx][ii][0] - system["orig"][0],
+                    system["coords"][f_idx][ii][1] - system["orig"][1],
+                    system["coords"][f_idx][ii][2] - system["orig"][2],
+                    system["spins"][f_idx][ii][0] / spins_norm[ii],
+                    system["spins"][f_idx][ii][1] / spins_norm[ii],
+                    system["spins"][f_idx][ii][2] / spins_norm[ii],
+                    spins_norm[ii],
+                )  # noqa: UP031
+            else:
+                ret += coord_fmt % (
+                    ii + 1,
+                    system["atom_types"][ii] + 1,
+                    system["coords"][f_idx][ii][0] - system["orig"][0],
+                    system["coords"][f_idx][ii][1] - system["orig"][1],
+                    system["coords"][f_idx][ii][2] - system["orig"][2],
+                    system["spins"][f_idx][ii][0],
+                    system["spins"][f_idx][ii][1],
+                    system["spins"][f_idx][ii][2] + 1,
+                    spins_norm[ii],
+                )  # noqa: UP031
+        else:
+            ret += coord_fmt % (
+                ii + 1,
+                system["atom_types"][ii] + 1,
+                system["coords"][f_idx][ii][0] - system["orig"][0],
+                system["coords"][f_idx][ii][1] - system["orig"][1],
+                system["coords"][f_idx][ii][2] - system["orig"][2],
+            )  # noqa: UP031
     return ret
 
 
diff --git a/dpdata/plugins/abacus.py b/dpdata/plugins/abacus.py
index cbd1475f2..b40d99668 100644
--- a/dpdata/plugins/abacus.py
+++ b/dpdata/plugins/abacus.py
@@ -20,7 +20,9 @@
 @Format.register("stru")
 class AbacusSTRUFormat(Format):
     def from_system(self, file_name, **kwargs):
-        return dpdata.abacus.scf.get_frame_from_stru(file_name)
+        data = dpdata.abacus.scf.get_frame_from_stru(file_name)
+        register_mag_data(data)
+        return data
 
     def to_system(self, data, file_name: FileType, frame_idx=0, **kwargs):
         """Dump the system into ABACUS STRU format file.
@@ -55,18 +57,32 @@ def register_mag_data(data):
             required=False,
             deepmd_name="spin",
         )
+        dpdata.System.register_data_type(dt)
         dpdata.LabeledSystem.register_data_type(dt)
-    if "mag_forces" in data:
+    if "force_mags" in data:
         dt = DataType(
-            "mag_forces",
+            "force_mags",
             np.ndarray,
             (Axis.NFRAMES, Axis.NATOMS, 3),
             required=False,
             deepmd_name="force_mag",
         )
+        dpdata.System.register_data_type(dt)
         dpdata.LabeledSystem.register_data_type(dt)
 
 
+def register_move_data(data):
+    if "move" in data:
+        dt = DataType(
+            "move",
+            np.ndarray,
+            (Axis.NFRAMES, Axis.NATOMS, 3),
+            required=False,
+            deepmd_name="move",
+        )
+        dpdata.System.register_data_type(dt)
+
+
 @Format.register("abacus/scf")
 @Format.register("abacus/pw/scf")
 @Format.register("abacus/lcao/scf")
@@ -75,6 +91,7 @@ class AbacusSCFFormat(Format):
     def from_labeled_system(self, file_name, **kwargs):
         data = dpdata.abacus.scf.get_frame(file_name)
         register_mag_data(data)
+        register_move_data(data)
         return data
 
 
@@ -86,6 +103,7 @@ class AbacusMDFormat(Format):
     def from_labeled_system(self, file_name, **kwargs):
         data = dpdata.abacus.md.get_frame(file_name)
         register_mag_data(data)
+        register_move_data(data)
         return data
 
 
@@ -97,4 +115,5 @@ class AbacusRelaxFormat(Format):
     def from_labeled_system(self, file_name, **kwargs):
         data = dpdata.abacus.relax.get_frame(file_name)
         register_mag_data(data)
+        register_move_data(data)
         return data
diff --git a/dpdata/plugins/amber.py b/dpdata/plugins/amber.py
index 361e0d8a7..c51af3465 100644
--- a/dpdata/plugins/amber.py
+++ b/dpdata/plugins/amber.py
@@ -136,8 +136,8 @@ def label(self, data: dict) -> dict:
         labeled_system = dpdata.LabeledSystem()
         with tempfile.TemporaryDirectory() as d:
             for ii, ss in enumerate(ori_system):
-                inp_fn = os.path.join(d, "%d.in" % ii)
-                out_fn = os.path.join(d, "%d.out" % ii)
+                inp_fn = os.path.join(d, "%d.in" % ii)  # noqa: UP031
+                out_fn = os.path.join(d, "%d.out" % ii)  # noqa: UP031
                 ss.to("sqm/in", inp_fn, **self.kwargs)
                 try:
                     sp.check_output(
diff --git a/dpdata/plugins/ase.py b/dpdata/plugins/ase.py
index 83a6bcf95..212d4af9b 100644
--- a/dpdata/plugins/ase.py
+++ b/dpdata/plugins/ase.py
@@ -175,7 +175,9 @@ def to_labeled_system(self, data, *args, **kwargs) -> list[ase.Atoms]:
                 cell=data["cells"][ii],
             )
 
-            results = {"energy": data["energies"][ii], "forces": data["forces"][ii]}
+            results = {"energy": data["energies"][ii]}
+            if "forces" in data:
+                results["forces"] = data["forces"][ii]
             if "virials" in data:
                 # convert to GPa as this is ase convention
                 # v_pref = 1 * 1e4 / 1.602176621e6
@@ -296,7 +298,10 @@ def from_labeled_system(
             dict_frames["energies"] = np.append(
                 dict_frames["energies"], tmp["energies"][0]
             )
-            dict_frames["forces"] = np.append(dict_frames["forces"], tmp["forces"][0])
+            if "forces" in tmp.keys() and "forces" in dict_frames.keys():
+                dict_frames["forces"] = np.append(
+                    dict_frames["forces"], tmp["forces"][0]
+                )
             if "virials" in tmp.keys() and "virials" in dict_frames.keys():
                 dict_frames["virials"] = np.append(
                     dict_frames["virials"], tmp["virials"][0]
@@ -305,7 +310,8 @@ def from_labeled_system(
         ## Correct the shape of numpy arrays
         dict_frames["cells"] = dict_frames["cells"].reshape(-1, 3, 3)
         dict_frames["coords"] = dict_frames["coords"].reshape(len(sub_traj), -1, 3)
-        dict_frames["forces"] = dict_frames["forces"].reshape(len(sub_traj), -1, 3)
+        if "forces" in dict_frames.keys():
+            dict_frames["forces"] = dict_frames["forces"].reshape(len(sub_traj), -1, 3)
         if "virials" in dict_frames.keys():
             dict_frames["virials"] = dict_frames["virials"].reshape(-1, 3, 3)
 
diff --git a/dpdata/plugins/deepmd.py b/dpdata/plugins/deepmd.py
index 1ed79f72e..2726e1d46 100644
--- a/dpdata/plugins/deepmd.py
+++ b/dpdata/plugins/deepmd.py
@@ -10,6 +10,7 @@
 import dpdata.deepmd.hdf5
 import dpdata.deepmd.mixed
 import dpdata.deepmd.raw
+from dpdata.data_type import Axis, DataType
 from dpdata.driver import Driver
 from dpdata.format import Format
 
@@ -17,10 +18,33 @@
     import h5py
 
 
+def register_spin():
+    dt = DataType(
+        "spins",
+        np.ndarray,
+        (Axis.NFRAMES, Axis.NATOMS, 3),
+        required=False,
+        deepmd_name="spin",
+    )
+    dpdata.System.register_data_type(dt)
+    dpdata.LabeledSystem.register_data_type(dt)
+
+    dt = DataType(
+        "force_mags",
+        np.ndarray,
+        (Axis.NFRAMES, Axis.NATOMS, 3),
+        required=False,
+        deepmd_name="force_mag",
+    )
+    dpdata.System.register_data_type(dt)
+    dpdata.LabeledSystem.register_data_type(dt)
+
+
 @Format.register("deepmd")
 @Format.register("deepmd/raw")
 class DeePMDRawFormat(Format):
     def from_system(self, file_name, type_map=None, **kwargs):
+        register_spin()
         return dpdata.deepmd.raw.to_system_data(
             file_name, type_map=type_map, labels=False
         )
@@ -30,6 +54,7 @@ def to_system(self, data, file_name, **kwargs):
         dpdata.deepmd.raw.dump(file_name, data)
 
     def from_labeled_system(self, file_name, type_map=None, **kwargs):
+        register_spin()
         return dpdata.deepmd.raw.to_system_data(
             file_name, type_map=type_map, labels=True
         )
@@ -41,6 +66,7 @@ def from_labeled_system(self, file_name, type_map=None, **kwargs):
 @Format.register("deepmd/comp")
 class DeePMDCompFormat(Format):
     def from_system(self, file_name, type_map=None, **kwargs):
+        register_spin()
         return dpdata.deepmd.comp.to_system_data(
             file_name, type_map=type_map, labels=False
         )
@@ -69,6 +95,7 @@ def to_system(self, data, file_name, set_size=5000, prec=np.float64, **kwargs):
         dpdata.deepmd.comp.dump(file_name, data, set_size=set_size, comp_prec=prec)
 
     def from_labeled_system(self, file_name, type_map=None, **kwargs):
+        register_spin()
         return dpdata.deepmd.comp.to_system_data(
             file_name, type_map=type_map, labels=True
         )
@@ -149,6 +176,7 @@ def mix_system(self, *system, type_map, **kwargs):
         return dpdata.deepmd.mixed.mix_system(*system, type_map=type_map, **kwargs)
 
     def from_multi_systems(self, directory, **kwargs):
+        register_spin()
         sys_dir = []
         for root, dirs, files in os.walk(directory):
             if (
@@ -204,6 +232,8 @@ def _from_system(
         """
         import h5py
 
+        register_spin()
+
         if isinstance(file_name, (h5py.Group, h5py.File)):
             return dpdata.deepmd.hdf5.to_system_data(
                 file_name, "", type_map=type_map, labels=labels
diff --git a/dpdata/plugins/gaussian.py b/dpdata/plugins/gaussian.py
index 55bee5a4c..9cba45989 100644
--- a/dpdata/plugins/gaussian.py
+++ b/dpdata/plugins/gaussian.py
@@ -109,8 +109,8 @@ def label(self, data: dict) -> dict:
         labeled_system = dpdata.LabeledSystem()
         with tempfile.TemporaryDirectory() as d:
             for ii, ss in enumerate(ori_system):
-                inp_fn = os.path.join(d, "%d.gjf" % ii)
-                out_fn = os.path.join(d, "%d.log" % ii)
+                inp_fn = os.path.join(d, "%d.gjf" % ii)  # noqa: UP031
+                out_fn = os.path.join(d, "%d.log" % ii)  # noqa: UP031
                 ss.to("gaussian/gjf", inp_fn, **self.kwargs)
                 try:
                     sp.check_output([*self.gaussian_exec.split(), inp_fn])
diff --git a/dpdata/plugins/lammps.py b/dpdata/plugins/lammps.py
index 0327d444b..c7e5c7653 100644
--- a/dpdata/plugins/lammps.py
+++ b/dpdata/plugins/lammps.py
@@ -2,8 +2,11 @@
 
 from typing import TYPE_CHECKING
 
+import numpy as np
+
 import dpdata.lammps.dump
 import dpdata.lammps.lmp
+from dpdata.data_type import Axis, DataType
 from dpdata.format import Format
 from dpdata.utils import open_file
 
@@ -11,6 +14,18 @@
     from dpdata.utils import FileType
 
 
+def register_spin(data):
+    if "spins" in data:
+        dt = DataType(
+            "spins",
+            np.ndarray,
+            (Axis.NFRAMES, Axis.NATOMS, 3),
+            required=False,
+            deepmd_name="spin",
+        )
+        dpdata.System.register_data_type(dt)
+
+
 @Format.register("lmp")
 @Format.register("lammps/lmp")
 class LAMMPSLmpFormat(Format):
@@ -18,7 +33,9 @@ class LAMMPSLmpFormat(Format):
     def from_system(self, file_name: FileType, type_map=None, **kwargs):
         with open_file(file_name) as fp:
             lines = [line.rstrip("\n") for line in fp]
-        return dpdata.lammps.lmp.to_system_data(lines, type_map)
+        data = dpdata.lammps.lmp.to_system_data(lines, type_map)
+        register_spin(data)
+        return data
 
     def to_system(self, data, file_name: FileType, frame_idx=0, **kwargs):
         """Dump the system in lammps data format.
@@ -45,7 +62,40 @@ def to_system(self, data, file_name: FileType, frame_idx=0, **kwargs):
 class LAMMPSDumpFormat(Format):
     @Format.post("shift_orig_zero")
     def from_system(
-        self, file_name, type_map=None, begin=0, step=1, unwrap=False, **kwargs
+        self,
+        file_name: str,
+        type_map: list[str] = None,
+        begin: int = 0,
+        step: int = 1,
+        unwrap: bool = False,
+        input_file: str = None,
+        **kwargs,
     ):
+        """Read the data from a lammps dump file.
+
+        Parameters
+        ----------
+        file_name : str
+            The dump file name
+        type_map : List[str], optional
+            The atom type list
+        begin : int, optional
+            The begin step
+        step : int, optional
+            The step
+        unwrap : bool, optional
+            Whether to unwrap the coordinates
+        input_file : str, optional
+            The input file name
+
+        Returns
+        -------
+        dict
+            The system data
+        """
         lines = dpdata.lammps.dump.load_file(file_name, begin=begin, step=step)
-        return dpdata.lammps.dump.system_data(lines, type_map, unwrap=unwrap)
+        data = dpdata.lammps.dump.system_data(
+            lines, type_map, unwrap=unwrap, input_file=input_file
+        )
+        register_spin(data)
+        return data
diff --git a/dpdata/plugins/n2p2.py b/dpdata/plugins/n2p2.py
index 28942ff5f..f9974fbab 100644
--- a/dpdata/plugins/n2p2.py
+++ b/dpdata/plugins/n2p2.py
@@ -85,18 +85,18 @@ def from_labeled_system(self, file_name: FileType, **kwargs):
                     energy = None
                 elif line.lower() == "end":
                     # If we are at the end of a section, process the section
-                    assert (
-                        len(coord) == len(atype) == len(force)
-                    ), "Number of atoms, atom types, and forces must match."
+                    assert len(coord) == len(atype) == len(force), (
+                        "Number of atoms, atom types, and forces must match."
+                    )
 
                     # Check if the number of atoms is consistent across all frames
                     natom = len(coord)
                     if natom0 is None:
                         natom0 = natom
                     else:
-                        assert (
-                            natom == natom0
-                        ), "The number of atoms in all frames must be the same."
+                        assert natom == natom0, (
+                            "The number of atoms in all frames must be the same."
+                        )
 
                     # Check if the number of atoms of each type is consistent across all frames
                     atype = np.array(atype)
@@ -108,9 +108,9 @@ def from_labeled_system(self, file_name: FileType, **kwargs):
                     if natoms0 is None:
                         natoms0 = natoms
                     else:
-                        assert (
-                            natoms == natoms0
-                        ), "The number of atoms of each type in all frames must be the same."
+                        assert natoms == natoms0, (
+                            "The number of atoms of each type in all frames must be the same."
+                        )
                     if atom_types0 is None:
                         atom_types0 = atype
                     atom_order = match_indices(atom_types0, atype)
diff --git a/dpdata/plugins/pwmat.py b/dpdata/plugins/pwmat.py
index ba3dab160..38a5bb297 100644
--- a/dpdata/plugins/pwmat.py
+++ b/dpdata/plugins/pwmat.py
@@ -31,11 +31,13 @@ def from_labeled_system(
             data["cells"],
             data["coords"],
             data["energies"],
-            data["forces"],
+            tmp_force,
             tmp_virial,
         ) = dpdata.pwmat.movement.get_frames(
             file_name, begin=begin, step=step, convergence_check=convergence_check
         )
+        if tmp_force is not None:
+            data["forces"] = tmp_force
         if tmp_virial is not None:
             data["virials"] = tmp_virial
         # scale virial to the unit of eV
diff --git a/dpdata/plugins/pymatgen.py b/dpdata/plugins/pymatgen.py
index 322298c3c..b8099a3ab 100644
--- a/dpdata/plugins/pymatgen.py
+++ b/dpdata/plugins/pymatgen.py
@@ -30,16 +30,15 @@ def to_system(self, data, **kwargs):
         """Convert System to Pymatgen Structure obj."""
         structures = []
         try:
-            from pymatgen.core import Structure
+            from pymatgen.core import Lattice, Structure
         except ModuleNotFoundError as e:
             raise ImportError("No module pymatgen.Structure") from e
 
-        species = []
-        for name, numb in zip(data["atom_names"], data["atom_numbs"]):
-            species.extend([name] * numb)
+        species = [data["atom_names"][tt] for tt in data["atom_types"]]
+        pbc = not (data.get("nopbc", False))
         for ii in range(data["coords"].shape[0]):
             structure = Structure(
-                data["cells"][ii],
+                Lattice(data["cells"][ii], pbc=[pbc] * 3),
                 species,
                 data["coords"][ii],
                 coords_are_cartesian=True,
diff --git a/dpdata/plugins/vasp.py b/dpdata/plugins/vasp.py
index 0160bde29..ca86ee193 100644
--- a/dpdata/plugins/vasp.py
+++ b/dpdata/plugins/vasp.py
@@ -7,6 +7,7 @@
 import dpdata.vasp.outcar
 import dpdata.vasp.poscar
 import dpdata.vasp.xml
+from dpdata.data_type import Axis, DataType
 from dpdata.format import Format
 from dpdata.utils import open_file, uniq_atom_names
 
@@ -14,6 +15,18 @@
     from dpdata.utils import FileType
 
 
+def register_move_data(data):
+    if "move" in data:
+        dt = DataType(
+            "move",
+            np.ndarray,
+            (Axis.NFRAMES, Axis.NATOMS, 3),
+            required=False,
+            deepmd_name="move",
+        )
+        dpdata.System.register_data_type(dt)
+
+
 @Format.register("poscar")
 @Format.register("contcar")
 @Format.register("vasp/poscar")
@@ -25,6 +38,7 @@ def from_system(self, file_name: FileType, **kwargs):
             lines = [line.rstrip("\n") for line in fp]
         data = dpdata.vasp.poscar.to_system_data(lines)
         data = uniq_atom_names(data)
+        register_move_data(data)
         return data
 
     def to_system(self, data, file_name: FileType, frame_idx=0, **kwargs):
@@ -81,7 +95,7 @@ def from_labeled_system(
             data["cells"],
             data["coords"],
             data["energies"],
-            data["forces"],
+            tmp_force,
             tmp_virial,
         ) = dpdata.vasp.outcar.get_frames(
             file_name,
@@ -90,6 +104,8 @@ def from_labeled_system(
             ml=ml,
             convergence_check=convergence_check,
         )
+        if tmp_force is not None:
+            data["forces"] = tmp_force
         if tmp_virial is not None:
             data["virials"] = tmp_virial
         # scale virial to the unit of eV
@@ -99,6 +115,7 @@ def from_labeled_system(
                 vol = np.linalg.det(np.reshape(data["cells"][ii], [3, 3]))
                 data["virials"][ii] *= v_pref * vol
         data = uniq_atom_names(data)
+        register_move_data(data)
         return data
 
 
@@ -135,4 +152,5 @@ def from_labeled_system(self, file_name, begin=0, step=1, **kwargs):
                 vol = np.linalg.det(np.reshape(data["cells"][ii], [3, 3]))
                 data["virials"][ii] *= v_pref * vol
         data = uniq_atom_names(data)
+        register_move_data(data)
         return data
diff --git a/dpdata/pwmat/atomconfig.py b/dpdata/pwmat/atomconfig.py
index 62eff77ca..5f01c8409 100644
--- a/dpdata/pwmat/atomconfig.py
+++ b/dpdata/pwmat/atomconfig.py
@@ -54,7 +54,7 @@ def to_system_data(lines):
 def from_system_data(system, f_idx=0, skip_zeros=True):
     ret = ""
     natoms = sum(system["atom_numbs"])
-    ret += "%d" % natoms
+    ret += "%d" % natoms  # noqa: UP031
     ret += "\n"
     ret += "LATTICE"
     ret += "\n"
@@ -83,7 +83,7 @@ def from_system_data(system, f_idx=0, skip_zeros=True):
     posi_list = []
     for jj, ii in zip(atomic_numbers, posis):
         ii = np.matmul(ii, np.linalg.inv(system["cells"][0]))
-        posi_list.append("%d %15.10f %15.10f %15.10f 1 1 1" % (jj, ii[0], ii[1], ii[2]))
+        posi_list.append("%d %15.10f %15.10f %15.10f 1 1 1" % (jj, ii[0], ii[1], ii[2]))  # noqa: UP031
     for kk in range(len(posi_list)):
         min = kk
         for jj in range(kk, len(posi_list)):
diff --git a/dpdata/pymatgen/structure.py b/dpdata/pymatgen/structure.py
index 36e411c02..1f74dbdd0 100644
--- a/dpdata/pymatgen/structure.py
+++ b/dpdata/pymatgen/structure.py
@@ -4,12 +4,19 @@
 
 
 def from_system_data(structure) -> dict:
-    symbols = [site.species_string for site in structure]
+    """Convert one pymatgen structure to dpdata's datadict."""
+    symbols = [ii.specie.symbol for ii in structure]
     atom_names = list(structure.symbol_set)
     atom_numbs = [symbols.count(symbol) for symbol in atom_names]
     atom_types = np.array([atom_names.index(symbol) for symbol in symbols]).astype(int)
     coords = structure.cart_coords
     cells = structure.lattice.matrix
+    if all(structure.pbc):
+        pbc = True
+    elif not any(structure.pbc):
+        pbc = False
+    else:
+        raise ValueError(f"Partial pbc condition {structure.pbc} is not supported")
 
     info_dict = {
         "atom_names": atom_names,
@@ -17,5 +24,7 @@ def from_system_data(structure) -> dict:
         "atom_types": atom_types,
         "coords": np.array([coords]),
         "cells": np.array([cells]),
+        "orig": np.zeros(3),
+        "nopbc": not pbc,
     }
     return info_dict
diff --git a/dpdata/qe/traj.py b/dpdata/qe/traj.py
index b4be303ad..4cb095500 100644
--- a/dpdata/qe/traj.py
+++ b/dpdata/qe/traj.py
@@ -78,7 +78,7 @@ def load_atom_names(lines, ntypes):
 def load_celldm(lines):
     celldm = np.zeros(6)
     for ii in range(6):
-        key = "celldm(%d)" % (ii + 1)
+        key = "celldm(%d)" % (ii + 1)  # noqa: UP031
         val = load_key(lines, key)
         if val is not None:
             celldm[ii] = float(val)
@@ -169,7 +169,7 @@ def load_energy(fname, begin=0, step=1):
     data = np.loadtxt(fname)
     steps = []
     for ii in data[begin::step, 0]:
-        steps.append("%d" % ii)
+        steps.append("%d" % ii)  # noqa: UP031
     with open_file(fname) as fp:
         while True:
             line = fp.readline()
diff --git a/dpdata/rdkit/sanitize.py b/dpdata/rdkit/sanitize.py
index 2b0d76634..1477918d2 100644
--- a/dpdata/rdkit/sanitize.py
+++ b/dpdata/rdkit/sanitize.py
@@ -367,10 +367,10 @@ def sanitize_mol(mol, verbose=False):
 # copy from FEprep
 def mol_edit_log(mol, i, j):
     if not mol.HasProp("edit"):
-        mol.SetProp("edit", "%d_%d" % (i, j))
+        mol.SetProp("edit", "%d_%d" % (i, j))  # noqa: UP031
     else:
         edited = mol.GetProp("edit")
-        mol.SetProp("edit", edited + ",%d_%d" % (i, j))
+        mol.SetProp("edit", edited + ",%d_%d" % (i, j))  # noqa: UP031
 
 
 def kekulize_aromatic_heterocycles(mol_in, assign_formal_charge=True, sanitize=True):
diff --git a/dpdata/stat.py b/dpdata/stat.py
index 5ec395708..0b04d5700 100644
--- a/dpdata/stat.py
+++ b/dpdata/stat.py
@@ -2,13 +2,14 @@
 
 from abc import ABCMeta, abstractmethod
 from functools import lru_cache
+from typing import Any
 
 import numpy as np
 
 from dpdata.system import LabeledSystem, MultiSystems
 
 
-def mae(errors: np.ndarray) -> np.float64:
+def mae(errors: np.ndarray) -> np.floating[Any]:
     """Compute the mean absolute error (MAE).
 
     Parameters
@@ -18,13 +19,13 @@ def mae(errors: np.ndarray) -> np.float64:
 
     Returns
     -------
-    np.float64
+    floating[Any]
         mean absolute error (MAE)
     """
     return np.mean(np.abs(errors))
 
 
-def rmse(errors: np.ndarray) -> np.float64:
+def rmse(errors: np.ndarray) -> np.floating[Any]:
     """Compute the root mean squared error (RMSE).
 
     Parameters
@@ -34,7 +35,7 @@ def rmse(errors: np.ndarray) -> np.float64:
 
     Returns
     -------
-    np.float64
+    floating[Any]
         root mean squared error (RMSE)
     """
     return np.sqrt(np.mean(np.square(errors)))
@@ -54,12 +55,12 @@ class ErrorsBase(metaclass=ABCMeta):
     SYSTEM_TYPE = object
 
     def __init__(self, system_1: SYSTEM_TYPE, system_2: SYSTEM_TYPE) -> None:
-        assert isinstance(
-            system_1, self.SYSTEM_TYPE
-        ), f"system_1 should be {self.SYSTEM_TYPE.__name__}"
-        assert isinstance(
-            system_2, self.SYSTEM_TYPE
-        ), f"system_2 should be {self.SYSTEM_TYPE.__name__}"
+        assert isinstance(system_1, self.SYSTEM_TYPE), (
+            f"system_1 should be {self.SYSTEM_TYPE.__name__}"
+        )
+        assert isinstance(system_2, self.SYSTEM_TYPE), (
+            f"system_2 should be {self.SYSTEM_TYPE.__name__}"
+        )
         self.system_1 = system_1
         self.system_2 = system_2
 
@@ -74,22 +75,22 @@ def f_errors(self) -> np.ndarray:
         """Force errors."""
 
     @property
-    def e_mae(self) -> np.float64:
+    def e_mae(self) -> np.floating[Any]:
         """Energy MAE."""
         return mae(self.e_errors)
 
     @property
-    def e_rmse(self) -> np.float64:
+    def e_rmse(self) -> np.floating[Any]:
         """Energy RMSE."""
         return rmse(self.e_errors)
 
     @property
-    def f_mae(self) -> np.float64:
+    def f_mae(self) -> np.floating[Any]:
         """Force MAE."""
         return mae(self.f_errors)
 
     @property
-    def f_rmse(self) -> np.float64:
+    def f_rmse(self) -> np.floating[Any]:
         """Force RMSE."""
         return rmse(self.f_errors)
 
@@ -115,7 +116,7 @@ class Errors(ErrorsBase):
     SYSTEM_TYPE = LabeledSystem
 
     @property
-    @lru_cache()
+    @lru_cache
     def e_errors(self) -> np.ndarray:
         """Energy errors."""
         assert isinstance(self.system_1, self.SYSTEM_TYPE)
@@ -123,7 +124,7 @@ def e_errors(self) -> np.ndarray:
         return self.system_1["energies"] - self.system_2["energies"]
 
     @property
-    @lru_cache()
+    @lru_cache
     def f_errors(self) -> np.ndarray:
         """Force errors."""
         assert isinstance(self.system_1, self.SYSTEM_TYPE)
@@ -152,7 +153,7 @@ class MultiErrors(ErrorsBase):
     SYSTEM_TYPE = MultiSystems
 
     @property
-    @lru_cache()
+    @lru_cache
     def e_errors(self) -> np.ndarray:
         """Energy errors."""
         assert isinstance(self.system_1, self.SYSTEM_TYPE)
@@ -165,7 +166,7 @@ def e_errors(self) -> np.ndarray:
         return np.concatenate(errors)
 
     @property
-    @lru_cache()
+    @lru_cache
     def f_errors(self) -> np.ndarray:
         """Force errors."""
         assert isinstance(self.system_1, self.SYSTEM_TYPE)
diff --git a/dpdata/system.py b/dpdata/system.py
index 08136cb91..cfc0f184e 100644
--- a/dpdata/system.py
+++ b/dpdata/system.py
@@ -5,21 +5,16 @@
 import hashlib
 import numbers
 import os
-import sys
 import warnings
 from copy import deepcopy
 from typing import (
     TYPE_CHECKING,
     Any,
     Iterable,
+    Literal,
     overload,
 )
 
-if sys.version_info >= (3, 8):
-    from typing import Literal
-else:
-    from typing_extensions import Literal
-
 import numpy as np
 
 import dpdata
@@ -226,7 +221,7 @@ def check_data(self):
             dd.check(self)
         if sum(self.get_atom_numbs()) != self.get_natoms():
             raise DataError(
-                "Sum of atom_numbs (%d) is not equal to natoms (%d)."
+                "Sum of atom_numbs (%d) is not equal to natoms (%d)."  # noqa: UP031
                 % (sum(self.get_atom_numbs()), self.get_natoms())
             )
 
@@ -281,8 +276,8 @@ def __str__(self):
         ret = "Data Summary"
         ret += "\nUnlabeled System"
         ret += "\n-------------------"
-        ret += "\nFrame Numbers     : %d" % self.get_nframes()
-        ret += "\nAtom Numbers      : %d" % self.get_natoms()
+        ret += "\nFrame Numbers     : %d" % self.get_nframes()  # noqa: UP031
+        ret += "\nAtom Numbers      : %d" % self.get_natoms()  # noqa: UP031
         ret += "\nElement List      :"
         ret += "\n-------------------"
         ret += "\n" + "  ".join(map(str, self.get_atom_names()))
@@ -1049,6 +1044,7 @@ def remove_atom_names(self, atom_names: str | list[str]):
             atom_idx = self.data["atom_types"] == idx
             removed_atom_idx.append(atom_idx)
         picked_atom_idx = ~np.any(removed_atom_idx, axis=0)
+        assert not isinstance(picked_atom_idx, np.bool_)
         new_sys = self.pick_atom_idx(picked_atom_idx)
         # let's remove atom_names
         # firstly, rearrange atom_names and put these atom_names in the end
@@ -1208,7 +1204,11 @@ class LabeledSystem(System):
     DTYPES: tuple[DataType, ...] = System.DTYPES + (
         DataType("energies", np.ndarray, (Axis.NFRAMES,), deepmd_name="energy"),
         DataType(
-            "forces", np.ndarray, (Axis.NFRAMES, Axis.NATOMS, 3), deepmd_name="force"
+            "forces",
+            np.ndarray,
+            (Axis.NFRAMES, Axis.NATOMS, 3),
+            required=False,
+            deepmd_name="force",
         ),
         DataType(
             "virials",
@@ -1243,8 +1243,8 @@ def __str__(self):
         ret = "Data Summary"
         ret += "\nLabeled System"
         ret += "\n-------------------"
-        ret += "\nFrame Numbers      : %d" % self.get_nframes()
-        ret += "\nAtom Numbers       : %d" % self.get_natoms()
+        ret += "\nFrame Numbers      : %d" % self.get_nframes()  # noqa: UP031
+        ret += "\nAtom Numbers       : %d" % self.get_natoms()  # noqa: UP031
         status = "Yes" if self.has_virial() else "No"
         ret += f"\nIncluding Virials  : {status}"
         ret += "\nElement List       :"
@@ -1268,13 +1268,17 @@ def __add__(self, others):
             raise RuntimeError("Unspported data structure")
         return self.__class__.from_dict({"data": self_copy.data})
 
+    def has_forces(self) -> bool:
+        return "forces" in self.data
+
     def has_virial(self) -> bool:
         # return ('virials' in self.data) and (len(self.data['virials']) > 0)
         return "virials" in self.data
 
     def affine_map_fv(self, trans, f_idx: int | numbers.Integral):
         assert np.linalg.det(trans) != 0
-        self.data["forces"][f_idx] = np.matmul(self.data["forces"][f_idx], trans)
+        if self.has_forces():
+            self.data["forces"][f_idx] = np.matmul(self.data["forces"][f_idx], trans)
         if self.has_virial():
             self.data["virials"][f_idx] = np.matmul(
                 trans.T, np.matmul(self.data["virials"][f_idx], trans)
@@ -1307,7 +1311,8 @@ def correction(self, hl_sys: LabeledSystem) -> LabeledSystem:
             raise RuntimeError("high_sys should be LabeledSystem")
         corrected_sys = self.copy()
         corrected_sys.data["energies"] = hl_sys.data["energies"] - self.data["energies"]
-        corrected_sys.data["forces"] = hl_sys.data["forces"] - self.data["forces"]
+        if "forces" in self.data and "forces" in hl_sys.data:
+            corrected_sys.data["forces"] = hl_sys.data["forces"] - self.data["forces"]
         if "virials" in self.data and "virials" in hl_sys.data:
             corrected_sys.data["virials"] = (
                 hl_sys.data["virials"] - self.data["virials"]
diff --git a/dpdata/unit.py b/dpdata/unit.py
index 09981b969..9c96827cd 100644
--- a/dpdata/unit.py
+++ b/dpdata/unit.py
@@ -2,13 +2,31 @@
 
 from abc import ABC
 
-from scipy import constants  # noqa: TID253
+# scipy physical constant version 2018
+physical_constants = {
+    "Avogadro constant": 6.02214076e23,  # mol^-1
+    "elementary charge": 1.602176634e-19,  # C
+    "atomic unit of length": 5.29177210903e-11,  # m
+    "atomic unit of energy": 4.3597447222071e-18,  # J
+    "Rydberg constant": 10973731.568160,  # m^-1
+    "Planck constant": 6.62607015e-34,  # J·s
+    "speed of light in vacuum": 299792458,  # m·s^-1
+}
+
+
+def scipy_constant_value(key: str) -> float:
+    return physical_constants[key]
+
 
-AVOGADRO = constants.Avogadro  # Avagadro constant
-ELE_CHG = constants.elementary_charge  # Elementary Charge, in C
-BOHR = constants.value("atomic unit of length")  # Bohr, in m
-HARTREE = constants.value("atomic unit of energy")  # Hartree, in Jole
-RYDBERG = constants.Rydberg * constants.h * constants.c  # Rydberg, in Jole
+AVOGADRO = scipy_constant_value("Avogadro constant")  # Avagadro constant
+ELE_CHG = scipy_constant_value("elementary charge")  # Elementary Charge, in C
+BOHR = scipy_constant_value("atomic unit of length")  # Bohr, in m
+HARTREE = scipy_constant_value("atomic unit of energy")  # Hartree, in Jole
+RYDBERG = (
+    scipy_constant_value("Rydberg constant")
+    * scipy_constant_value("Planck constant")
+    * scipy_constant_value("speed of light in vacuum")
+)  # Rydberg, in Jole
 
 # energy conversions
 econvs = {
diff --git a/dpdata/utils.py b/dpdata/utils.py
index 8942bd54d..58a908cc7 100644
--- a/dpdata/utils.py
+++ b/dpdata/utils.py
@@ -2,14 +2,9 @@
 
 import io
 import os
-import sys
 from contextlib import contextmanager
-from typing import TYPE_CHECKING, Generator, overload
+from typing import TYPE_CHECKING, Generator, Literal, overload
 
-if sys.version_info >= (3, 8):
-    from typing import Literal
-else:
-    from typing_extensions import Literal
 import numpy as np
 
 from dpdata.periodic_table import Element
diff --git a/dpdata/vasp/outcar.py b/dpdata/vasp/outcar.py
index 0fa4cb68e..08b8fce96 100644
--- a/dpdata/vasp/outcar.py
+++ b/dpdata/vasp/outcar.py
@@ -161,9 +161,9 @@ def analyze_block(lines, ntot, nelm, ml=False):
                 not lines[idx + in_kB_index].split()[0:2] == ["in", "kB"]
             ):
                 in_kB_index += 1
-            assert idx + in_kB_index < len(
-                lines
-            ), 'ERROR: "in kB" is not found in OUTCAR. Unable to extract virial.'
+            assert idx + in_kB_index < len(lines), (
+                'ERROR: "in kB" is not found in OUTCAR. Unable to extract virial.'
+            )
             tmp_v = [float(ss) for ss in lines[idx + in_kB_index].split()[2:8]]
             virial = np.zeros([3, 3])
             virial[0][0] = tmp_v[0]
diff --git a/dpdata/vasp/poscar.py b/dpdata/vasp/poscar.py
index 102e79041..78b8dbbeb 100644
--- a/dpdata/vasp/poscar.py
+++ b/dpdata/vasp/poscar.py
@@ -4,13 +4,22 @@
 import numpy as np
 
 
-def _to_system_data_lower(lines, cartesian=True):
+def _to_system_data_lower(lines, cartesian=True, selective_dynamics=False):
+    def move_flag_mapper(flag):
+        if flag == "T":
+            return True
+        elif flag == "F":
+            return False
+        else:
+            raise RuntimeError(f"Invalid move flag: {flag}")
+
     """Treat as cartesian poscar."""
     system = {}
     system["atom_names"] = [str(ii) for ii in lines[5].split()]
     system["atom_numbs"] = [int(ii) for ii in lines[6].split()]
     scale = float(lines[1])
     cell = []
+    move_flags = []
     for ii in range(2, 5):
         boxv = [float(jj) for jj in lines[ii].split()]
         boxv = np.array(boxv) * scale
@@ -19,12 +28,21 @@ def _to_system_data_lower(lines, cartesian=True):
     natoms = sum(system["atom_numbs"])
     coord = []
     for ii in range(8, 8 + natoms):
-        tmpv = [float(jj) for jj in lines[ii].split()[:3]]
+        tmp = lines[ii].split()
+        tmpv = [float(jj) for jj in tmp[:3]]
         if cartesian:
             tmpv = np.array(tmpv) * scale
         else:
             tmpv = np.matmul(np.array(tmpv), system["cells"][0])
         coord.append(tmpv)
+        if selective_dynamics:
+            if len(tmp) == 6:
+                move_flags.append(list(map(move_flag_mapper, tmp[3:])))
+            else:
+                raise RuntimeError(
+                    f"Invalid move flags, should be 6 columns, got {tmp}"
+                )
+
     system["coords"] = [np.array(coord)]
     system["orig"] = np.zeros(3)
     atom_types = []
@@ -34,12 +52,18 @@ def _to_system_data_lower(lines, cartesian=True):
     system["atom_types"] = np.array(atom_types, dtype=int)
     system["cells"] = np.array(system["cells"])
     system["coords"] = np.array(system["coords"])
+    if move_flags:
+        move_flags = np.array(move_flags, dtype=bool)
+        move_flags = move_flags.reshape((1, natoms, 3))
+        system["move"] = np.array(move_flags, dtype=bool)
     return system
 
 
 def to_system_data(lines):
     # remove the line that has 'selective dynamics'
+    selective_dynamics = False
     if lines[7][0] == "S" or lines[7][0] == "s":
+        selective_dynamics = True
         lines.pop(7)
     is_cartesian = lines[7][0] in ["C", "c", "K", "k"]
     if not is_cartesian:
@@ -47,7 +71,7 @@ def to_system_data(lines):
             raise RuntimeError(
                 "seem not to be a valid POSCAR of vasp 5.x, may be a POSCAR of vasp 4.x?"
             )
-    return _to_system_data_lower(lines, is_cartesian)
+    return _to_system_data_lower(lines, is_cartesian, selective_dynamics)
 
 
 def from_system_data(system, f_idx=0, skip_zeros=True):
@@ -55,7 +79,7 @@ def from_system_data(system, f_idx=0, skip_zeros=True):
     for ii, name in zip(system["atom_numbs"], system["atom_names"]):
         if ii == 0:
             continue
-        ret += "%s%d " % (name, ii)
+        ret += "%s%d " % (name, ii)  # noqa: UP031
     ret += "\n"
     ret += "1.0\n"
     for ii in system["cells"][f_idx]:
@@ -70,8 +94,12 @@ def from_system_data(system, f_idx=0, skip_zeros=True):
     for ii in system["atom_numbs"]:
         if ii == 0:
             continue
-        ret += "%d " % ii
+        ret += "%d " % ii  # noqa: UP031
     ret += "\n"
+    move = system.get("move", None)
+    if move is not None and len(move) > 0:
+        ret += "Selective Dynamics\n"
+
     # should use Cartesian for VESTA software
     ret += "Cartesian\n"
     atype = system["atom_types"]
@@ -81,9 +109,26 @@ def from_system_data(system, f_idx=0, skip_zeros=True):
     sort_idx = np.lexsort((np.arange(len(atype)), atype))
     atype = atype[sort_idx]
     posis = posis[sort_idx]
+    if move is not None and len(move) > 0:
+        move = move[f_idx][sort_idx]
+
+    if isinstance(move, np.ndarray):
+        move = move.tolist()
+
     posi_list = []
-    for ii in posis:
-        posi_list.append(f"{ii[0]:15.10f} {ii[1]:15.10f} {ii[2]:15.10f}")
+    for idx in range(len(posis)):
+        ii_posi = posis[idx]
+        line = f"{ii_posi[0]:15.10f} {ii_posi[1]:15.10f} {ii_posi[2]:15.10f}"
+        if move is not None and len(move) > 0:
+            move_flags = move[idx]
+            if not isinstance(move_flags, list) or len(move_flags) != 3:
+                raise RuntimeError(
+                    f"Invalid move flags: {move_flags}, should be a list of 3 bools"
+                )
+            line += " " + " ".join("T" if flag else "F" for flag in move_flags)
+
+        posi_list.append(line)
+
     posi_list.append("")
     ret += "\n".join(posi_list)
     return ret
diff --git a/dpdata/vasp/xml.py b/dpdata/vasp/xml.py
index 352b107ed..91c29bf9c 100755
--- a/dpdata/vasp/xml.py
+++ b/dpdata/vasp/xml.py
@@ -7,9 +7,9 @@
 
 
 def check_name(item, name):
-    assert (
-        item.attrib["name"] == name
-    ), "item attrib '{}' dose not math required '{}'".format(item.attrib["name"], name)
+    assert item.attrib["name"] == name, (
+        "item attrib '{}' dose not math required '{}'".format(item.attrib["name"], name)
+    )
 
 
 def get_varray(varray):
@@ -58,7 +58,7 @@ def formulate_config(eles, types, posi, cell, ener, forc, strs_):
     natoms = len(types)
     ntypes = len(eles)
     ret = ""
-    ret += "#N %d %d\n" % (natoms, ntypes - 1)
+    ret += "#N %d %d\n" % (natoms, ntypes - 1)  # noqa: UP031
     ret += "#C "
     for ii in eles:
         ret += " " + ii
@@ -73,7 +73,7 @@ def formulate_config(eles, types, posi, cell, ener, forc, strs_):
     ret += "#F\n"
     for ii in range(natoms):
         sp = np.matmul(cell.T, posi[ii])
-        ret += "%d" % (types[ii] - 1)
+        ret += "%d" % (types[ii] - 1)  # noqa: UP031
         ret += f" {sp[0]:12.6f} {sp[1]:12.6f} {sp[2]:12.6f}"
         ret += f" {forc[ii][0]:12.6f} {forc[ii][1]:12.6f} {forc[ii][2]:12.6f}"
         ret += "\n"
diff --git a/dpdata/xyz/quip_gap_xyz.py b/dpdata/xyz/quip_gap_xyz.py
index b23b27e07..cba971e47 100644
--- a/dpdata/xyz/quip_gap_xyz.py
+++ b/dpdata/xyz/quip_gap_xyz.py
@@ -161,7 +161,7 @@ def handle_single_xyz_frame(lines):
                         list(filter(bool, field_dict["virial"].split(" ")))
                     ).reshape(3, 3)
                 ]
-            ).astype("float32")
+            ).astype(np.float64)
         else:
             virials = None
 
@@ -175,10 +175,10 @@ def handle_single_xyz_frame(lines):
                     3, 3
                 )
             ]
-        ).astype("float32")
-        info_dict["coords"] = np.array([coords_array]).astype("float32")
-        info_dict["energies"] = np.array([field_dict["energy"]]).astype("float32")
-        info_dict["forces"] = np.array([force_array]).astype("float32")
+        ).astype(np.float64)
+        info_dict["coords"] = np.array([coords_array]).astype(np.float64)
+        info_dict["energies"] = np.array([field_dict["energy"]]).astype(np.float64)
+        info_dict["forces"] = np.array([force_array]).astype(np.float64)
         if virials is not None:
             info_dict["virials"] = virials
         info_dict["orig"] = np.zeros(3)
diff --git a/pyproject.toml b/pyproject.toml
index 4bf2e64b7..52c47804e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,6 @@ authors = [
 ]
 license = {file = "LICENSE"}
 classifiers = [
-    "Programming Language :: Python :: 3.7",
     "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
@@ -28,7 +27,7 @@ dependencies = [
     'importlib_metadata>=1.4; python_version < "3.8"',
     'typing_extensions; python_version < "3.8"',
 ]
-requires-python = ">=3.7"
+requires-python = ">=3.8"
 readme = "README.md"
 keywords = ["lammps", "vasp", "deepmd-kit"]
 
@@ -37,14 +36,11 @@ Homepage = "https://github.com/deepmodeling/dpdata"
 documentation = "https://docs.deepmodeling.com/projects/dpdata"
 repository = "https://github.com/deepmodeling/dpdata"
 
-[project.entry-points.console_scripts]
+[project.scripts]
 dpdata = "dpdata.cli:dpdata_cli"
 
 [project.optional-dependencies]
 test = [
-    # https://github.com/materialsproject/pymatgen/issues/3882
-    # https://github.com/kuelumbus/rdkit-pypi/issues/102
-    "numpy<2",
 ]
 ase = ['ase']
 amber = [
@@ -55,10 +51,10 @@ pymatgen = ['pymatgen']
 docs = [
     'sphinx',
     'recommonmark',
-    'sphinx_rtd_theme>=1.0.0rc1',
+    'sphinx-book-theme',
     'numpydoc',
     'myst-parser',
-    'deepmodeling-sphinx>=0.1.1',
+    'deepmodeling-sphinx>=0.3.0',
     'sphinx-argparse<0.5.0',
     'rdkit',
     'jupyterlite-sphinx',
diff --git a/tests/abacus.scf/STRU.ch4 b/tests/abacus.scf/STRU.ch4
index cf97747aa..bc33cbe94 100644
--- a/tests/abacus.scf/STRU.ch4
+++ b/tests/abacus.scf/STRU.ch4
@@ -18,11 +18,11 @@ Cartesian 				#Cartesian(Unit is LATTICE_CONSTANT)
 C 						#Name of element	
 0.0						#Magnetic for this element.
 1				#Number of atoms
-0.981274803    0.861285385    0.838442496 0 0 0
+0.981274803    0.861285385    0.838442496 1 1 1
 H
 0.0
 4
 1.023557202    0.758025625    0.66351336 0 0 0
-0.78075702    0.889445935    0.837363468 0 0 0
-1.064091613    1.043438905    0.840995502 0 0 0
-1.039321214    0.756530859    1.009609207 0 0 0
+0.78075702    0.889445935    0.837363468 1 0 1
+1.064091613    1.043438905    0.840995502 1 0 1
+1.039321214    0.756530859    1.009609207 0 1 1
diff --git a/tests/abacus.scf/stru_test b/tests/abacus.scf/stru_test
index 22d619c93..c3a1917db 100644
--- a/tests/abacus.scf/stru_test
+++ b/tests/abacus.scf/stru_test
@@ -22,11 +22,11 @@ Cartesian    # Cartesian(Unit is LATTICE_CONSTANT)
 C
 0.0
 1
-5.192682633809 4.557725978258 4.436846615358 1 1 1
+5.192682633809 4.557725978258 4.436846615358 1 1 1 mag 0.000000000000 0.000000000000 0.000000000000
 H
 0.0
 4
-5.416431453540 4.011298860305 3.511161492417 1 1 1
-4.131588222365 4.706745191323 4.431136645083 1 1 1
-5.630930319126 5.521640894956 4.450356541303 1 1 1
-5.499851012568 4.003388899277 5.342621842622 1 1 1
+5.416431453540 4.011298860305 3.511161492417 0 0 0 mag 0.000000000000 0.000000000000 0.000000000000
+4.131588222365 4.706745191323 4.431136645083 1 0 1 mag 0.000000000000 0.000000000000 0.000000000000
+5.630930319126 5.521640894956 4.450356541303 1 0 1 mag 0.000000000000 0.000000000000 0.000000000000
+5.499851012568 4.003388899277 5.342621842622 0 1 1 mag 0.000000000000 0.000000000000 0.000000000000
diff --git a/tests/abacus.spin/INPUT b/tests/abacus.spin/INPUT.scf.nspin2
similarity index 92%
rename from tests/abacus.spin/INPUT
rename to tests/abacus.spin/INPUT.scf.nspin2
index bc0b1df0d..0bdd93039 100644
--- a/tests/abacus.spin/INPUT
+++ b/tests/abacus.spin/INPUT.scf.nspin2
@@ -1,5 +1,5 @@
 INPUT_PARAMETERS
-suffix	ABACUS
+suffix	ABACUS-nspin2
 calculation	scf
 ecutwfc	100
 scf_thr	1e-07
@@ -14,7 +14,7 @@ basis_type	lcao
 symmetry	0
 noncolin	1
 lspinorb	0
-nspin	4
+nspin	2
 out_mul	true
 sc_mag_switch	1
 decay_grad_switch	1
diff --git a/tests/abacus.spin/OUT.ABACUS-nspin2/running_scf.log b/tests/abacus.spin/OUT.ABACUS-nspin2/running_scf.log
new file mode 100644
index 000000000..60cdf1563
--- /dev/null
+++ b/tests/abacus.spin/OUT.ABACUS-nspin2/running_scf.log
@@ -0,0 +1,565 @@
+                                                                                     
+                              ABACUS v3.7.1
+
+               Atomic-orbital Based Ab-initio Computation at UStc                    
+
+                     Website: http://abacus.ustc.edu.cn/                             
+               Documentation: https://abacus.deepmodeling.com/                       
+                  Repository: https://github.com/abacusmodeling/abacus-develop       
+                              https://github.com/deepmodeling/abacus-develop         
+                      Commit: 1a7a3158b (Fri Aug 23 00:52:25 2024 +0800)
+
+    Start Time is Fri Aug 23 16:06:02 2024
+                                                                                     
+ ------------------------------------------------------------------------------------
+
+ READING GENERAL INFORMATION
+                           global_out_dir = OUT.ABACUS/
+                           global_in_card = INPUT
+                               pseudo_dir = ./
+                              orbital_dir = ./
+                                    DRANK = 1
+                                    DSIZE = 16
+                                   DCOLOR = 1
+                                    GRANK = 1
+                                    GSIZE = 1
+
+
+
+
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ |                                                                    |
+ | Reading atom information in unitcell:                              |
+ | From the input file and the structure file we know the number of   |
+ | different elments in this unitcell, then we list the detail        |
+ | information for each element, especially the zeta and polar atomic |
+ | orbital number for each element. The total atom number is counted. |
+ | We calculate the nearest atom distance for each atom and show the  |
+ | Cartesian and Direct coordinates for each atom. We list the file   |
+ | address for atomic orbitals. The volume and the lattice vectors    |
+ | in real and reciprocal space is also shown.                        |
+ |                                                                    |
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+
+ READING UNITCELL INFORMATION
+                                    ntype = 1
+                  lattice constant (Bohr) = 1.88028
+              lattice constant (Angstrom) = 0.995
+
+ READING ATOM TYPE 1
+                               atom label = Fe
+                      L=0, number of zeta = 4
+                      L=1, number of zeta = 2
+                      L=2, number of zeta = 2
+                      L=3, number of zeta = 1
+             number of atom for this type = 2
+               magnetization of element 1 = [ 0, 0, 2.4 ]
+
+                        TOTAL ATOM NUMBER = 2
+DIRECT COORDINATES
+    atom                   x                   y                   z     mag                  vx                  vy                  vz
+taud_Fe1            0.0500000000        0.1000000000        0.1500000000  2.4000        0.0000000000        0.0000000000        0.0000000000
+taud_Fe2            0.5000000000        0.5000000000        0.5000000000  2.4000        0.0000000000        0.0000000000        0.0000000000
+
+
+
+                          Volume (Bohr^3) = 150.259
+                             Volume (A^3) = 22.266
+
+ Lattice vectors: (Cartesian coordinate: in unit of a_0)
+             +2.82743                  +0                  +0
+                   +0            +2.82743                  +0
+                   +0                  +0            +2.82743
+ Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)
+            +0.353679                  -0                  +0
+                   -0           +0.353679                  -0
+                   +0                  -0           +0.353679
+ The esolver type has been set to : ksdft_pw
+
+ RUNNING WITH DEVICE  : CPU / Intel(R) Xeon(R) Platinum 8255C CPU @ 2.50GHz
+
+
+
+
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ |                                                                    |
+ | Reading pseudopotentials files:                                    |
+ | The pseudopotential file is in UPF format. The 'NC' indicates that |
+ | the type of pseudopotential is 'norm conserving'. Functional of    |
+ | exchange and correlation is decided by 4 given parameters in UPF   |
+ | file.  We also read in the 'core correction' if there exists.      |
+ | Also we can read the valence electrons number and the maximal      |
+ | angular momentum used in this pseudopotential. We also read in the |
+ | trail wave function, trail atomic density and local-pseudopotential|
+ | on logrithmic grid. The non-local pseudopotential projector is also|
+ | read in if there is any.                                           |
+ |                                                                    |
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+
+                PAO radial cut off (Bohr) = 15
+
+ Read in pseudopotential file is Fe.upf
+                     pseudopotential type = NC
+          exchange-correlation functional = PBE
+                 nonlocal core correction = 1
+                        valence electrons = 16
+                                     lmax = 2
+                           number of zeta = 6
+                     number of projectors = 10
+                           L of projector = 0
+                           L of projector = 0
+                           L of projector = 1
+                           L of projector = 1
+                           L of projector = 1
+                           L of projector = 1
+                           L of projector = 2
+                           L of projector = 2
+                           L of projector = 2
+                           L of projector = 2
+     initial pseudo atomic orbital number = 40
+                                   NLOCAL = 108
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ Warning: the number of valence electrons in pseudopotential > 8 for Fe: [Ar] 3d6 4s2
+ Pseudopotentials with additional electrons can yield (more) accurate outcomes, but may be less efficient.
+ If you're confident that your chosen pseudopotential is appropriate, you can safely ignore this warning.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+                  
+
+
+
+
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ |                                                                    |
+ | Setup plane waves of charge/potential:                             |
+ | Use the energy cutoff and the lattice vectors to generate the      |
+ | dimensions of FFT grid. The number of FFT grid on each processor   |
+ | is 'nrxx'. The number of plane wave basis in reciprocal space is   |
+ | different for charege/potential and wave functions. We also set    |
+ | the 'sticks' for the parallel of FFT. The number of plane waves    |
+ | is 'npw' in each processor.                                        |
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+
+
+ SETUP THE PLANE WAVE BASIS
+ energy cutoff for charge/potential (unit:Ry) = 320
+            fft grid for charge/potential = [ 32, 32, 32 ]
+                        fft grid division = [ 1, 1, 1 ]
+        big fft grid for charge/potential = [ 32, 32, 32 ]
+                                     nbxx = 2048
+                                     nrxx = 2048
+
+ SETUP PLANE WAVES FOR CHARGE/POTENTIAL
+                    number of plane waves = 14531
+                         number of sticks = 725
+
+ PARALLEL PW FOR CHARGE/POTENTIAL
+     PROC   COLUMNS(POT)             PW
+        1             45            909
+        2             45            909
+        3             46            908
+        4             46            908
+        5             45            909
+        6             45            909
+        7             45            909
+        8             45            909
+        9             45            909
+       10             46            908
+       11             46            908
+       12             46            908
+       13             45            907
+       14             45            907
+       15             45            907
+       16             45            907
+ --------------- sum -------------------
+       16            725          14531
+                            number of |g| = 165
+                                  max |g| = 28.6453
+                                  min |g| = 0.500354
+
+----------- Double Check Mixing Parameters Begin ------------
+mixing_type: broyden
+mixing_beta: 0.4
+mixing_gg0: 1
+mixing_gg0_min: 0.1
+mixing_beta_mag: 1.6
+mixing_gg0_mag: 0
+mixing_ndim: 10
+----------- Double Check Mixing Parameters End ------------
+
+ SETUP THE ELECTRONS NUMBER
+            electron number of element Fe = 16
+      total electron number of element Fe = 32
+            AUTOSET number of electrons:  = 32
+ DONE : SETUP UNITCELL Time : 0.316237 (SEC)
+
+
+
+
+
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ |                                                                    |
+ | Setup K-points                                                     |
+ | We setup the k-points according to input parameters.               |
+ | The reduced k-points are set according to symmetry operations.     |
+ | We treat the spin as another set of k-points.                      |
+ |                                                                    |
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+
+
+ SETUP K-POINTS
+                                    nspin = 4
+                   Input type of k points = Monkhorst-Pack(Gamma)
+                                   nkstot = 216
+K-POINTS DIRECT COORDINATES
+ KPOINTS    DIRECT_X    DIRECT_Y    DIRECT_Z  WEIGHT
+
+ DONE : INIT K-POINTS Time : 0.376998 (SEC)
+
+
+
+
+
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ |                                                                    |
+ | Setup plane waves of wave functions:                               |
+ | Use the energy cutoff and the lattice vectors to generate the      |
+ | dimensions of FFT grid. The number of FFT grid on each processor   |
+ | is 'nrxx'. The number of plane wave basis in reciprocal space is   |
+ | different for charege/potential and wave functions. We also set    |
+ | the 'sticks' for the parallel of FFT. The number of plane wave of  |
+ | each k-point is 'npwk[ik]' in each processor                       |
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+
+
+ SETUP PLANE WAVES FOR WAVE FUNCTIONS
+     energy cutoff for wavefunc (unit:Ry) = 80
+              fft grid for wave functions = [ 32, 32, 32 ]
+                    number of plane waves = 3071
+                         number of sticks = 253
+
+ PARALLEL PW FOR WAVE FUNCTIONS
+     PROC   COLUMNS(POT)             PW
+        1             16            192
+        2             16            192
+        3             16            192
+        4             16            192
+        5             15            191
+        6             15            191
+        7             16            194
+        8             16            192
+        9             16            192
+       10             16            192
+       11             16            192
+       12             16            192
+       13             16            192
+       14             16            192
+       15             16            192
+       16             15            191
+ --------------- sum -------------------
+       16            253           3071
+ DONE : INIT PLANEWAVE Time : 0.380345 (SEC)
+
+                           occupied bands = 32
+                                   NBANDS = 52
+
+ SETUP NONLOCAL PSEUDOPOTENTIALS IN PLANE WAVE BASIS
+ Fe non-local projectors:
+ projector 1 L=0
+ projector 2 L=0
+ projector 3 L=1
+ projector 4 L=1
+ projector 5 L=1
+ projector 6 L=1
+ projector 7 L=2
+ projector 8 L=2
+ projector 9 L=2
+ projector 10 L=2
+      TOTAL NUMBER OF NONLOCAL PROJECTORS = 68
+ DONE : LOCAL POTENTIAL Time : 0.40963 (SEC)
+
+
+ Init Non-Local PseudoPotential table : 
+ Init Non-Local-Pseudopotential done.
+ DONE : NON-LOCAL POTENTIAL Time : 0.443773 (SEC)
+
+                                     npwx = 128
+
+ Warning_Memory_Consuming allocated:  Psi_PW 43.875 MB
+
+ Make real space PAO into reciprocal space.
+       max mesh points in Pseudopotential = 1425
+     dq(describe PAO in reciprocal space) = 0.01
+                                    max q = 1078
+
+ number of pseudo atomic orbitals for Fe is 6
+ DONE : INIT BASIS Time : 0.620862 (SEC)
+
+
+ -------------------------------------------
+ SELF-CONSISTENT
+ -------------------------------------------
+                                 init_chg = atomic
+ DONE : INIT SCF Time : 0.761759 (SEC)
+
+
+ PW ALGORITHM --------------- ION=   1  ELEC=   1--------------------------------
+total magnetism (Bohr mag/cell)	-8.56154e-05	-3.16838e-05	3.87672
+       absolute magnetism (Bohr mag/cell) = 3.90219
+ Notice: Threshold on eigenvalues was too large.
+ hsover_error=0.32 > DRHO=0.0810864
+ Origin diag_ethr = 0.01
+ New    diag_ethr = 0.000253395
+total magnetism (Bohr mag/cell)	-9.09654e-05	-3.00029e-05	3.87658
+       absolute magnetism (Bohr mag/cell) = 3.90286
+
+ Density error is 0.0818013187196
+                          Error Threshold = 0.000253394954616
+----------------------------------------------------------
+     Energy           Rydberg                 eV          
+----------------------------------------------------------
+ E_KohnSham     -501.1638820263      -6818.6844273569     
+ E_Harris       -501.1650591987      -6818.7004436094     
+ E_Fermi        1.4203911993         19.3254136999        
+----------------------------------------------------------
+
+
+ PW ALGORITHM --------------- ION=   1  ELEC=   2--------------------------------
+total magnetism (Bohr mag/cell)	-2.92622585913e-05	-0.000266812474323	3.65800660509
+       absolute magnetism (Bohr mag/cell) = 3.88387266047
+
+ Density error is 0.0165562104235
+                          Error Threshold = 0.000255629120999
+----------------------------------------------------------
+     Energy           Rydberg                 eV          
+----------------------------------------------------------
+ E_KohnSham     -501.1785311471      -6818.8837388713     
+ E_Harris       -501.4729927600      -6822.8900946492     
+ E_Fermi        1.4123735299         19.2163277107        
+----------------------------------------------------------
+
+
+ PW ALGORITHM --------------- ION=   1  ELEC=   3--------------------------------
+total magnetism (Bohr mag/cell)	-1.41288005953e-05	-0.000388528787066	3.72548472333
+       absolute magnetism (Bohr mag/cell) = 3.91023806289
+
+ Density error is 0.00318633115585
+                          Error Threshold = 5.17381575736e-05
+----------------------------------------------------------
+     Energy           Rydberg                 eV          
+----------------------------------------------------------
+ E_KohnSham     -501.1785047719      -6818.8833800180     
+ E_Harris       -501.0742053019      -6817.4643129270     
+ E_Fermi        1.4117399370         19.2077072372        
+----------------------------------------------------------
+
+
+ PW ALGORITHM --------------- ION=   1  ELEC=   4--------------------------------
+total magnetism (Bohr mag/cell)	2.26526732328e-05	-0.000450436470697	3.73775099313
+       absolute magnetism (Bohr mag/cell) = 3.92308219614
+
+ Density error is 0.000165938860835
+                          Error Threshold = 9.95728486202e-06
+----------------------------------------------------------
+     Energy           Rydberg                 eV          
+----------------------------------------------------------
+ E_KohnSham     -501.1804562269      -6818.9099309251     
+ E_Harris       -501.0912756169      -6817.6965664782     
+ E_Fermi        1.4131895976         19.2274308815        
+----------------------------------------------------------
+
+
+ PW ALGORITHM --------------- ION=   1  ELEC=   5--------------------------------
+
+total magnetism (Bohr mag/cell)	-0.000682480084963	-0.000407948602903	4.47635324671
+       absolute magnetism (Bohr mag/cell) = 4.82248432181
+
+ Density error is 6.61858307885e-09
+                          Error Threshold = 5.40730613116e-10
+----------------------------------------------------------
+     Energy           Rydberg                 eV          
+----------------------------------------------------------
+ E_KohnSham     -501.1664531630      -6818.7194094666     
+ E_KS(sigma->0) -501.1658847659      -6818.7116760274     
+ E_Harris       -501.1680065331      -6818.7405441510     
+ E_band         -30.1701507306       -410.4859594555      
+ E_one_elec     -214.9321986489      -2924.3025852934     
+ E_Hartree      128.2299640744       1744.6581657467      
+ E_xc           -71.3503189041       -970.7708912124      
+ E_Ewald        -343.2422281761      -4670.0500974110     
+ E_entropy(-TS) -0.0011367942        -0.0154668785        
+ E_descf        0.0000000000         0.0000000000         
+ E_exx          0.0000000000         0.0000000000         
+ E_Fermi        1.4108073605         19.1950188828        
+----------------------------------------------------------
+
+-------------------------------------------------------------------------------------------
+ Total Magnetism (uB)                                                                      
+-------------------------------------------------------------------------------------------
+                         Fe         2.4000001004 
+                         Fe         2.3999994597 
+-------------------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------------------
+ Magnetic force (eV/uB)                                                                    
+-------------------------------------------------------------------------------------------
+                         Fe         -0.3669618965 
+                         Fe         -0.3669821632 
+-------------------------------------------------------------------------------------------
+
+Charge and Mag of Atom 0
+ Orbital 1 Charge: 1.99985926204 Mag: -8.50835272933e-09 -8.84955876115e-09 2.84869655853e-06
+ Orbital 2 Charge: 5.99701367096 Mag: -2.18088393737e-07 -2.40618119919e-07 0.000105269162551
+ Orbital 3 Charge: 6.37379729477 Mag: -4.581340875e-08 7.66548679323e-08 2.39989198254
+Sum of atom 0 is: 14.3706702278 -2.72410155216e-07 -1.72812810748e-07 2.4000001004
+Charge and Mag of Atom 1
+ Orbital 1 Charge: 1.99985925962 Mag: -6.91586723071e-09 -8.93979594505e-09 2.85090960195e-06
+ Orbital 2 Charge: 5.99701367924 Mag: -2.16539560943e-07 -2.55979271358e-07 0.000105265428294
+ Orbital 3 Charge: 6.37381969542 Mag: -9.45758719447e-08 3.50118887944e-08 2.3998913434
+Sum of atom 1 is: 14.3706926343 -3.18031300118e-07 -2.29907178509e-07 2.39999945974
+
+ charge density convergence is achieved
+ final etot is -6818.7194095 eV
+ EFERMI = 19.195018883 eV
+
+------------------------------------------------------------------------------------------
+ TOTAL-FORCE (eV/Angstrom)                                                                
+------------------------------------------------------------------------------------------
+                       Fe1        -0.0000234480         0.0000848730         0.0003959639 
+                       Fe2         0.0000234480        -0.0000848730        -0.0003959639 
+------------------------------------------------------------------------------------------
+
+ --------------------------------------------
+ !FINAL_ETOT_IS -6818.719409466637 eV
+ --------------------------------------------
+
+
+TIME STATISTICS
+----------------------------------------------------------------------------
+    CLASS_NAME               NAME             TIME/s  CALLS   AVG/s  PER/%  
+----------------------------------------------------------------------------
+                  total                       563.56 15       37.57  100.00 
+ Driver           reading                     0.11   1        0.11   0.02   
+ Input_Conv       Convert                     0.00   1        0.00   0.00   
+ Driver           driver_line                 563.45 1        563.45 99.98  
+ UnitCell         check_tau                   0.00   1        0.00   0.00   
+ PW_Basis_Sup     setuptransform              0.06   1        0.06   0.01   
+ PW_Basis_Sup     distributeg                 0.00   1        0.00   0.00   
+ mymath           heapsort                    0.00   3        0.00   0.00   
+ PW_Basis_K       setuptransform              0.00   1        0.00   0.00   
+ PW_Basis_K       distributeg                 0.00   1        0.00   0.00   
+ PW_Basis         setup_struc_factor          0.00   1        0.00   0.00   
+ ppcell_vnl       init                        0.00   1        0.00   0.00   
+ ppcell_vl        init_vloc                   0.02   1        0.02   0.00   
+ ppcell_vnl       init_vnl                    0.03   1        0.03   0.01   
+ WF_atomic        init_at_1                   0.18   1        0.18   0.03   
+ wavefunc         wfcinit                     0.00   1        0.00   0.00   
+ Ions             opt_ions                    562.87 1        562.87 99.88  
+ ESolver_KS_PW    runner                      562.87 1        562.87 99.88  
+ H_Ewald_pw       compute_ewald               0.00   1        0.00   0.00   
+ Charge           set_rho_core                0.01   1        0.01   0.00   
+ PW_Basis_Sup     recip2real                  0.03   202      0.00   0.00   
+ PW_Basis_Sup     gathers_scatterp            0.01   202      0.00   0.00   
+ Charge           atomic_rho                  0.00   1        0.00   0.00   
+ Potential        init_pot                    0.01   1        0.01   0.00   
+ Potential        update_from_charge          0.07   13       0.01   0.01   
+ Potential        cal_fixed_v                 0.00   1        0.00   0.00   
+ PotLocal         cal_fixed_v                 0.00   1        0.00   0.00   
+ Potential        cal_v_eff                   0.06   13       0.00   0.01   
+ H_Hartree_pw     v_hartree                   0.00   13       0.00   0.00   
+ PW_Basis_Sup     real2recip                  0.02   247      0.00   0.00   
+ PW_Basis_Sup     gatherp_scatters            0.02   247      0.00   0.00   
+ PotXC            cal_v_eff                   0.06   13       0.00   0.01   
+ XC_Functional    v_xc                        0.06   13       0.00   0.01   
+ Potential        interpolate_vrs             0.00   13       0.00   0.00   
+ OnsiteProj       init_k_stage0               0.08   1        0.08   0.01   
+ Charge_Mixing    init_mixing                 0.00   2        0.00   0.00   
+ ESolver_KS_PW    hamilt2density              184.77 13       14.21  32.79  
+ HSolverPW        solve                       287.77 21       13.70  51.06  
+ pp_cell_vnl      getvnl                      4.21   21384    0.00   0.75   
+ Structure_Factor get_sk                      2.26   492264   0.00   0.40   
+ OnsiteProj       getvnl                      3.45   21384    0.00   0.61   
+ OnsiteProj       init_k                      4.11   26136    0.00   0.73   
+ OnsiteProj       init_k_stage1               1.73   26136    0.00   0.31   
+ OnsiteProj       init_k_stage2               2.37   26136    0.00   0.42   
+ WF_atomic        atomic_wfc                  0.10   216      0.00   0.02   
+ DiagoDavid       diag_mock                   259.79 4536     0.06   46.10  
+ DiagoDavid       first                       83.93  4536     0.02   14.89  
+ David            spsi_func                   1.10   1075598  0.00   0.20   
+ DiagoDavid       SchmidtOrth                 14.54  537799   0.00   2.58   
+ David            hpsi_func                   149.60 14979    0.01   26.55  
+ Operator         hPsi                        396.41 31827    0.01   70.34  
+ Operator         EkineticPW                  2.03   31827    0.00   0.36   
+ Operator         VeffPW                      344.42 31827    0.01   61.12  
+ PW_Basis_K       recip2real                  183.20 3119822  0.00   32.51  
+ PW_Basis_K       gathers_scatterp            160.47 3119822  0.00   28.47  
+ PW_Basis_K       real2recip                  155.75 2827790  0.00   27.64  
+ PW_Basis_K       gatherp_scatters            136.28 2827790  0.00   24.18  
+ Operator         NonlocalPW                  42.43  31827    0.00   7.53   
+ Nonlocal         add_nonlocal_pp             34.57  31827    0.00   6.13   
+ Operator         OnsiteProjPW                7.38   31827    0.00   1.31   
+ OnsiteProj       add_onsite_proj             7.35   31827    0.00   1.30   
+ OnsiteProj       overlap                     5.91   36579    0.00   1.05   
+ DiagoDavid       cal_elem                    35.38  14979    0.00   6.28   
+ DiagoDavid       diag_zhegvx                 51.18  14979    0.00   9.08   
+ DiagoDavid       cal_grad                    101.08 10443    0.01   17.94  
+ DiagoDavid       check_update                0.01   10443    0.00   0.00   
+ DiagoDavid       last                        1.75   4536     0.00   0.31   
+ ElecStatePW      psiToRho                    18.60  13       1.43   3.30   
+ Charge_Mixing    get_drho                    0.01   13       0.00   0.00   
+ Charge_Mixing    inner_product_recip_rho     0.00   13       0.00   0.00   
+ Charge           mix_rho                     0.01   10       0.00   0.00   
+ Charge           Broyden_mixing              0.01   10       0.00   0.00   
+ Charge_Mixing    inner_product_recip_hartree 0.00   48       0.00   0.00   
+ SpinConstrain    cal_mw_from_lambda          377.05 86       4.38   66.90  
+ DiagoIterAssist  diag_responce               266.63 16848    0.02   47.31  
+ DiagoIterAssist  diagH_LAPACK                12.16  16848    0.00   2.16   
+ OnsiteProj       cal_occupation              0.06   1        0.06   0.01   
+ ModuleIO         write_istate_info           0.07   1        0.07   0.01   
+----------------------------------------------------------------------------
+
+
+ NAME-------------------------|MEMORY(MB)--------
+                         total          786.3
+                        Psi_PW          668.8
+                    DAV::basis          12.38
+                     DAV::hpsi          12.38
+                     DAV::spsi          12.38
+                      DAV::hcc          10.56
+                      DAV::scc          10.56
+                      DAV::vcc          10.56
+                      MTransOp          10.56
+                  PW_B_K::gcar          9.645
+                   PW_B_K::gk2          3.215
+                 DiagSub::temp          3.096
+                      VNL::tab          2.632
+                      VNL::vkb          2.024
+            Nonlocal<PW>::becp          1.727
+              Nonlocal<PW>::ps          1.727
+                   VNL::tab_at          1.579
+                     FFT::grid          1.000
+                      Chg::rho          1.000
+                 Chg::rho_save          1.000
+                     Pot::veff          1.000
+              Pot::veff_smooth          1.000
+ -------------   < 1.0 MB has been ignored ----------------
+ ----------------------------------------------------------
+
+ Start  Time  : Fri Aug 23 16:06:02 2024
+ Finish Time  : Fri Aug 23 16:15:04 2024
+ Total  Time  : 0 h 9 mins 2 secs 
diff --git a/tests/abacus.spin/STRU.spin b/tests/abacus.spin/STRU.spin
new file mode 100644
index 000000000..340430c90
--- /dev/null
+++ b/tests/abacus.spin/STRU.spin
@@ -0,0 +1,24 @@
+ATOMIC_SPECIES
+Fe 55.845 Fe.upf
+
+NUMERICAL_ORBITAL
+Fe_gga_10au_200.0Ry_4s2p2d1f.orb
+
+LATTICE_CONSTANT
+1.880277359
+LATTICE_VECTORS
+    2.8274254848     0.0000000000     0.0000000000 #latvec1
+    0.0000000000     2.8274254848     0.0000000000 #latvec2
+    0.0000000000     0.0000000000     2.8274254848 #latvec3
+
+ATOMIC_POSITIONS
+Direct
+
+Fe #label
+1 #magnetism
+4 #number of atoms
+    0.0000000000     0.000000000     0.000000000 mag 0 0 2
+    0.1000000000     0.1000000000     0.1000000000 mag 3
+    0.2000000000     0.2000000000     0.2000000000 mag 3 angle1 90
+    0.3000000000     0.3000000000     0.3000000000 mag 3 4 0 angle1 90 angle2 90
+
diff --git a/tests/comp_sys.py b/tests/comp_sys.py
index a3a916a09..4a38e4d30 100644
--- a/tests/comp_sys.py
+++ b/tests/comp_sys.py
@@ -118,15 +118,15 @@ def _make_comp_ms_test_func(comp_sys_test_func):
     """
 
     def comp_ms_test_func(iobj):
-        assert hasattr(iobj, "ms_1") and hasattr(
-            iobj, "ms_2"
-        ), "Multi-system objects must be present"
+        assert hasattr(iobj, "ms_1") and hasattr(iobj, "ms_2"), (
+            "Multi-system objects must be present"
+        )
         iobj.assertEqual(len(iobj.ms_1), len(iobj.ms_2))
         keys = [ii.formula for ii in iobj.ms_1]
         keys_2 = [ii.formula for ii in iobj.ms_2]
-        assert sorted(keys) == sorted(
-            keys_2
-        ), f"Keys of two MS are not equal: {keys} != {keys_2}"
+        assert sorted(keys) == sorted(keys_2), (
+            f"Keys of two MS are not equal: {keys} != {keys_2}"
+        )
         for kk in keys:
             iobj.system_1 = iobj.ms_1[kk]
             iobj.system_2 = iobj.ms_2[kk]
@@ -197,17 +197,17 @@ def test_is_nopbc(self):
 
 class MSAllIsPBC:
     def test_is_pbc(self):
-        assert hasattr(self, "ms_1") and hasattr(
-            self, "ms_2"
-        ), "Multi-system objects must be present and iterable"
+        assert hasattr(self, "ms_1") and hasattr(self, "ms_2"), (
+            "Multi-system objects must be present and iterable"
+        )
         self.assertTrue(all([not ss.nopbc for ss in self.ms_1]))
         self.assertTrue(all([not ss.nopbc for ss in self.ms_2]))
 
 
 class MSAllIsNoPBC:
     def test_is_nopbc(self):
-        assert hasattr(self, "ms_1") and hasattr(
-            self, "ms_2"
-        ), "Multi-system objects must be present and iterable"
+        assert hasattr(self, "ms_1") and hasattr(self, "ms_2"), (
+            "Multi-system objects must be present and iterable"
+        )
         self.assertTrue(all([ss.nopbc for ss in self.ms_1]))
         self.assertTrue(all([ss.nopbc for ss in self.ms_2]))
diff --git a/tests/lammps/in.lmp b/tests/lammps/in.lmp
new file mode 100644
index 000000000..02f655fa5
--- /dev/null
+++ b/tests/lammps/in.lmp
@@ -0,0 +1,2 @@
+compute 	spin all property/atom sp spx spy spz fmx fmy fmz fx fy fz
+dump            dpgen_dump all custom 10 traj.dump id type x y z c_spin[1] c_spin[2] c_spin[3] c_spin[4] c_spin[5] c_spin[6] c_spin[7] c_spin[8] c_spin[9] c_spin[10]
\ No newline at end of file
diff --git a/tests/lammps/spin.lmp b/tests/lammps/spin.lmp
new file mode 100644
index 000000000..d662cc2cd
--- /dev/null
+++ b/tests/lammps/spin.lmp
@@ -0,0 +1,12 @@
+
+2 atoms
+2 atom types
+   0.0000000000    2.5243712000 xlo xhi
+   0.0000000000    2.0430257000 ylo yhi
+   0.0000000000    2.2254033000 zlo zhi
+   1.2621856000    1.2874292000    0.7485898000 xy xz yz
+
+Atoms # atomic
+
+     1      1    0.0000000000    0.0000000000    0.0000000000    0.6000000000    0.8000000000    0.0000000000    5.0000000000
+     2      2    1.2621856000    0.7018028000    0.5513885000    0.0000000000    0.8000000000    0.6000000000    5.0000000000
diff --git a/tests/lammps/traj.dump b/tests/lammps/traj.dump
new file mode 100644
index 000000000..e1682983a
--- /dev/null
+++ b/tests/lammps/traj.dump
@@ -0,0 +1,52 @@
+ITEM: TIMESTEP
+0
+ITEM: NUMBER OF ATOMS
+17
+ITEM: BOX BOUNDS xy xz yz pp pp pp
+-4.0080511965879438e-02 5.7039029418994556e+00 -5.9179115295410201e-03
+1.4436085788922526e-02 5.6674744441011660e+00 -1.1487414836883500e-02
+7.8239288740356017e-03 5.6734038274259646e+00 6.8277359008788905e-04
+ITEM: ATOMS id type x y z c_spin[1] c_spin[2] c_spin[3] c_spin[4] c_spin[5] c_spin[6] c_spin[7] c_spin[8] c_spin[9] c_spin[10]
+1 1 0.00141160 5.64868599 0.01005602 1.54706291 0.00000000 0.00000000 1.00000000 -1.40772100 -2.03739417 -1522.64797384 -0.00397809 -0.00190426 -0.00743976
+2 1 5.65283939 5.57449025 2.84281508 1.54412869 0.00000000 0.00000000 1.00000000 7.75304092 6.48949619 -1512.84926162 -0.00637234 -0.00733168 0.00661107
+3 1 0.00066480 2.78022036 0.01010716 1.54612979 0.00000000 0.00000000 1.00000000 -0.93618575 1.92206111 -1520.80305011 -0.00316673 0.00177893 -0.00744575
+4 1 5.65233666 2.85374747 2.84289453 1.54439093 0.00000000 0.00000000 1.00000000 8.11012818 -6.49922039 -1514.31557088 -0.00569217 0.00741000 0.00640353
+5 1 2.82063515 5.64869321 0.01007552 1.54714250 0.00000000 0.00000000 1.00000000 2.49070852 -2.14456666 -1523.53038650 0.00478410 -0.00213962 -0.00751154
+6 1 2.89579803 5.57439179 2.84287630 1.54415032 0.00000000 0.00000000 1.00000000 -8.03062338 6.63950296 -1513.41291897 0.00440396 -0.00717185 0.00633657
+7 1 2.82151287 2.78010538 0.01016303 1.54619615 0.00000000 0.00000000 1.00000000 2.71859584 1.98482729 -1521.34149633 0.00533453 0.00194532 -0.00745901
+8 1 2.89637049 2.85377083 2.84297332 1.54440023 0.00000000 0.00000000 1.00000000 -7.76758760 -6.67134514 -1514.43304618 0.00505040 0.00743195 0.00630302
+9 1 1.41106492 1.38817482 1.72302072 1.18134529 0.00000000 0.00000000 1.00000000 0.27170165 -0.00426695 -444.22843899 0.00100237 -0.00002725 -0.03385965
+10 1 1.41105247 1.38807861 3.96314606 1.18153407 0.00000000 0.00000000 1.00000000 -0.07722674 0.01368756 -337.08703133 -0.00066982 0.00007487 0.07887183
+11 1 1.41105864 4.21395432 1.43987180 1.71989299 0.00000000 0.00000000 1.00000000 -0.01511106 0.00320081 -1653.34500916 0.00010421 0.00007248 0.00634401
+12 1 1.41104843 4.21387554 4.24576823 1.71989825 0.00000000 0.00000000 1.00000000 -0.71645898 0.05923960 -1640.68070568 -0.00117959 0.00006676 -0.01467806
+13 1 4.27433865 1.38779084 1.43977211 1.72010048 0.00000000 0.00000000 1.00000000 0.45899480 0.03956420 -1653.36356942 0.00051885 0.00002313 0.00911600
+14 1 4.27436799 1.38772964 4.24586490 1.72010133 0.00000000 0.00000000 1.00000000 0.38385331 0.07301994 -1642.06086017 -0.00002034 0.00010335 -0.01688908
+15 1 4.27435427 4.21452597 1.39359689 1.65590121 0.00000000 0.00000000 1.00000000 -0.01658773 -0.06159007 -1659.12744163 0.00006470 -0.00006420 -0.01342935
+16 1 4.27434583 4.21455469 4.29208004 1.65592002 0.00000000 0.00000000 1.00000000 -0.15590720 -0.03252166 -1654.84697132 -0.00066755 -0.00003915 -0.00482188
+17 2 1.41096761 1.38958048 0.01029027 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00048351 -0.00022876 -0.00645195
+ITEM: TIMESTEP
+10
+ITEM: NUMBER OF ATOMS
+17
+ITEM: BOX BOUNDS xy xz yz pp pp pp
+-4.0080511965879438e-02 5.7039029418994556e+00 -5.9179115295410201e-03
+1.4436085788922526e-02 5.6674744441011660e+00 -1.1487414836883500e-02
+7.8239288740356017e-03 5.6734038274259646e+00 6.8277359008788905e-04
+ITEM: ATOMS id type x y z c_spin[1] c_spin[2] c_spin[3] c_spin[4] c_spin[5] c_spin[6] c_spin[7] c_spin[8] c_spin[9] c_spin[10]
+1 1 0.00037565 5.64900783 0.00994919 1.20356102 0.17466098 0.84115562 -0.51181127 -77.61471611 -389.41594243 234.29512368 0.00514290 -0.02481576 0.01063015
+2 1 5.65299480 5.57370279 2.84182058 1.17910451 0.85296110 0.48195380 -0.20044424 -311.75775120 -175.76677913 79.45225558 -0.01239308 -0.00914070 0.01933082
+3 1 0.00076668 2.78053566 0.01181481 1.20779106 0.33415542 0.49831517 -0.80001384 -163.88630094 -248.58823709 387.72415159 -0.01738465 -0.02878227 0.01503087
+4 1 5.65188602 2.85285383 2.84413423 1.20124335 -0.83536303 -0.20314716 0.51078356 399.86863784 90.34522869 -236.39221701 0.02327635 -0.00046572 -0.00138388
+5 1 2.82101290 5.64942265 0.01091135 1.34670883 -0.98528016 0.07078135 -0.15560530 902.73741755 -62.52279896 140.44423419 0.01500524 0.00581151 0.00468238
+6 1 2.89400594 5.57477971 2.84333235 1.25424131 -0.94587492 0.11487066 0.30352161 528.43507318 -60.32699018 -171.89948334 -0.00478280 0.00069273 -0.00496159
+7 1 2.82260306 2.78052696 0.00917962 1.17249564 -0.99589145 0.06282562 -0.06521619 374.56568243 -26.39431071 20.98877908 0.01464486 -0.01010131 -0.00993410
+8 1 2.89632273 2.85545549 2.84070353 1.24297017 -0.44008251 -0.42493729 0.79104721 240.05525392 236.02796206 -448.18443804 0.00137705 0.01258804 -0.01817420
+9 1 1.41117683 1.38867159 1.72266429 1.19059484 0.71251804 -0.69714805 -0.07938914 -309.93474514 293.96860716 19.98886311 -0.03871152 0.00854863 -0.02757569
+10 1 1.41176544 1.38873530 3.96470435 1.17564502 -0.51932602 -0.74875017 0.41191463 181.72443401 263.91689829 -132.94216896 0.00122847 0.01674701 0.02707109
+11 1 1.41085716 4.21342650 1.43850987 1.19874662 -0.51890828 0.82913822 0.20800000 237.52969259 -379.65100512 -93.16140268 0.01185986 -0.01872789 0.00032128
+12 1 1.41088045 4.21340876 4.24487134 1.20157661 -0.86390154 -0.04516556 -0.50163154 388.97171693 21.75492170 227.68580658 0.02074490 0.00756366 0.01937948
+13 1 4.27525485 1.38812593 1.43912039 1.23209806 0.55809649 0.81404794 0.16079259 -335.92026314 -484.87463129 -91.14464759 -0.03675759 -0.03549076 0.00310277
+14 1 4.27483864 1.38696457 4.24782541 1.18431742 0.00519166 -0.92210080 0.38691492 -4.73957478 407.09534135 -171.59043210 -0.00911750 0.04394272 -0.01683249
+15 1 4.27528588 4.21463764 1.39334117 1.17456490 -0.93713453 -0.09927163 0.33455046 397.32993706 40.92599847 -141.68618750 0.01918926 -0.00534149 -0.01906574
+16 1 4.27407834 4.21327842 4.29226033 1.31499905 -0.21350543 -0.97682201 -0.01530327 180.98908307 848.25344747 12.36402507 0.00492895 0.04383813 0.00955221
+17 2 1.40675897 1.38612182 0.01000617 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00174929 -0.00686653 -0.01117336
diff --git a/tests/lammps/traj_partial_spin.dump b/tests/lammps/traj_partial_spin.dump
new file mode 100644
index 000000000..5ee20855a
--- /dev/null
+++ b/tests/lammps/traj_partial_spin.dump
@@ -0,0 +1,52 @@
+ITEM: TIMESTEP
+0
+ITEM: NUMBER OF ATOMS
+17
+ITEM: BOX BOUNDS xy xz yz pp pp pp
+-4.0080511965879438e-02 5.7039029418994556e+00 -5.9179115295410201e-03
+1.4436085788922526e-02 5.6674744441011660e+00 -1.1487414836883500e-02
+7.8239288740356017e-03 5.6734038274259646e+00 6.8277359008788905e-04
+ITEM: ATOMS id type x y z c_spin[1] c_spin[2] c_spin[3] c_spin[4] c_spin[5] c_spin[6] c_spin[7] c_spin[8] c_spin[9] c_spin[10]
+1 1 0.00141160 5.64868599 0.01005602 1.54706291 0.00000000 0.00000000 1.00000000 -1.40772100 -2.03739417 -1522.64797384 -0.00397809 -0.00190426 -0.00743976
+2 1 5.65283939 5.57449025 2.84281508 1.54412869 0.00000000 0.00000000 1.00000000 7.75304092 6.48949619 -1512.84926162 -0.00637234 -0.00733168 0.00661107
+3 1 0.00066480 2.78022036 0.01010716 1.54612979 0.00000000 0.00000000 1.00000000 -0.93618575 1.92206111 -1520.80305011 -0.00316673 0.00177893 -0.00744575
+4 1 5.65233666 2.85374747 2.84289453 1.54439093 0.00000000 0.00000000 1.00000000 8.11012818 -6.49922039 -1514.31557088 -0.00569217 0.00741000 0.00640353
+5 1 2.82063515 5.64869321 0.01007552 1.54714250 0.00000000 0.00000000 1.00000000 2.49070852 -2.14456666 -1523.53038650 0.00478410 -0.00213962 -0.00751154
+6 1 2.89579803 5.57439179 2.84287630 1.54415032 0.00000000 0.00000000 1.00000000 -8.03062338 6.63950296 -1513.41291897 0.00440396 -0.00717185 0.00633657
+7 1 2.82151287 2.78010538 0.01016303 1.54619615 0.00000000 0.00000000 1.00000000 2.71859584 1.98482729 -1521.34149633 0.00533453 0.00194532 -0.00745901
+8 1 2.89637049 2.85377083 2.84297332 1.54440023 0.00000000 0.00000000 1.00000000 -7.76758760 -6.67134514 -1514.43304618 0.00505040 0.00743195 0.00630302
+9 1 1.41106492 1.38817482 1.72302072 1.18134529 0.00000000 0.00000000 1.00000000 0.27170165 -0.00426695 -444.22843899 0.00100237 -0.00002725 -0.03385965
+10 1 1.41105247 1.38807861 3.96314606 1.18153407 0.00000000 0.00000000 1.00000000 -0.07722674 0.01368756 -337.08703133 -0.00066982 0.00007487 0.07887183
+11 1 1.41105864 4.21395432 1.43987180 1.71989299 0.00000000 0.00000000 1.00000000 -0.01511106 0.00320081 -1653.34500916 0.00010421 0.00007248 0.00634401
+12 1 1.41104843 4.21387554 4.24576823 1.71989825 0.00000000 0.00000000 1.00000000 -0.71645898 0.05923960 -1640.68070568 -0.00117959 0.00006676 -0.01467806
+13 1 4.27433865 1.38779084 1.43977211 1.72010048 0.00000000 0.00000000 1.00000000 0.45899480 0.03956420 -1653.36356942 0.00051885 0.00002313 0.00911600
+14 1 4.27436799 1.38772964 4.24586490 1.72010133 0.00000000 0.00000000 1.00000000 0.38385331 0.07301994 -1642.06086017 -0.00002034 0.00010335 -0.01688908
+15 1 4.27435427 4.21452597 1.39359689 1.65590121 0.00000000 0.00000000 1.00000000 -0.01658773 -0.06159007 -1659.12744163 0.00006470 -0.00006420 -0.01342935
+16 1 4.27434583 4.21455469 4.29208004 1.65592002 0.00000000 0.00000000 1.00000000 -0.15590720 -0.03252166 -1654.84697132 -0.00066755 -0.00003915 -0.00482188
+17 2 1.41096761 1.38958048 0.01029027 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00048351 -0.00022876 -0.00645195
+ITEM: TIMESTEP
+10
+ITEM: NUMBER OF ATOMS
+17
+ITEM: BOX BOUNDS xy xz yz pp pp pp
+-4.0080511965879438e-02 5.7039029418994556e+00 -5.9179115295410201e-03
+1.4436085788922526e-02 5.6674744441011660e+00 -1.1487414836883500e-02
+7.8239288740356017e-03 5.6734038274259646e+00 6.8277359008788905e-04
+ITEM: ATOMS id type x y z 
+1 1 0.00037565 5.64900783 0.00994919 
+2 1 5.65299480 5.57370279 2.84182058 
+3 1 0.00076668 2.78053566 0.01181481 
+4 1 5.65188602 2.85285383 2.84413423 
+5 1 2.82101290 5.64942265 0.01091135 
+6 1 2.89400594 5.57477971 2.84333235 
+7 1 2.82260306 2.78052696 0.00917962 
+8 1 2.89632273 2.85545549 2.84070353 
+9 1 1.41117683 1.38867159 1.72266429 
+10 1 1.41176544 1.38873530 3.96470435
+11 1 1.41085716 4.21342650 1.43850987
+12 1 1.41088045 4.21340876 4.24487134
+13 1 4.27525485 1.38812593 1.43912039
+14 1 4.27483864 1.38696457 4.24782541
+15 1 4.27528588 4.21463764 1.39334117
+16 1 4.27407834 4.21327842 4.29226033
+17 2 1.40675897 1.38612182 0.01000617
diff --git a/tests/plugin/pyproject.toml b/tests/plugin/pyproject.toml
index 7ce1f854a..3e01f27cb 100644
--- a/tests/plugin/pyproject.toml
+++ b/tests/plugin/pyproject.toml
@@ -11,7 +11,7 @@ dependencies = [
     'dpdata',
 ]
 readme = "README.md"
-requires-python = ">=3.7"
+requires-python = ">=3.8"
 
 [project.entry-points.'dpdata.plugins']
 random = "dpdata_plugin_test:ep"
diff --git a/tests/poscars/POSCAR.oh.err1 b/tests/poscars/POSCAR.oh.err1
new file mode 100644
index 000000000..1b7521e6e
--- /dev/null
+++ b/tests/poscars/POSCAR.oh.err1
@@ -0,0 +1,11 @@
+Cubic BN
+   3.57
+ 0.00 0.50 0.50
+ 0.45 0.00 0.50
+ 0.55 0.51 0.00
+ O H
+   1 1
+Selective dynamics
+Cartesian
+ 0.00 0.00 0.00 T T F
+ 0.25 0.25 0.25 F F
diff --git a/tests/poscars/POSCAR.oh.err2 b/tests/poscars/POSCAR.oh.err2
new file mode 100644
index 000000000..ed52d4d08
--- /dev/null
+++ b/tests/poscars/POSCAR.oh.err2
@@ -0,0 +1,11 @@
+Cubic BN
+   3.57
+ 0.00 0.50 0.50
+ 0.45 0.00 0.50
+ 0.55 0.51 0.00
+ O H
+   1 1
+Selective dynamics
+Cartesian
+ 0.00 0.00 0.00 T T F
+ 0.25 0.25 0.25 a T F
diff --git a/tests/poscars/poscar_ref_oh.py b/tests/poscars/poscar_ref_oh.py
index 2d29aeeb6..39299ff88 100644
--- a/tests/poscars/poscar_ref_oh.py
+++ b/tests/poscars/poscar_ref_oh.py
@@ -32,7 +32,7 @@ def test_cell(self):
                     self.system.data["cells"][0][ii][jj],
                     ovito_cell[ii][jj],
                     places=6,
-                    msg="cell[%d][%d] failed" % (ii, jj),
+                    msg="cell[%d][%d] failed" % (ii, jj),  # noqa: UP031
                 )
 
     def test_frame(self):
@@ -48,5 +48,5 @@ def test_frame(self):
                     self.system.data["coords"][0][ii][jj],
                     ovito_posis[ii][jj],
                     places=6,
-                    msg="posis[%d][%d] failed" % (ii, jj),
+                    msg="posis[%d][%d] failed" % (ii, jj),  # noqa: UP031
                 )
diff --git a/tests/pwmat/config_ref_ch4.py b/tests/pwmat/config_ref_ch4.py
index 72499398e..bb76f6cff 100644
--- a/tests/pwmat/config_ref_ch4.py
+++ b/tests/pwmat/config_ref_ch4.py
@@ -35,7 +35,7 @@ def test_cell(self):
                     self.system.data["cells"][0][ii][jj],
                     ovito_cell[ii][jj],
                     places=6,
-                    msg="cell[%d][%d] failed" % (ii, jj),
+                    msg="cell[%d][%d] failed" % (ii, jj),  # noqa: UP031
                 )
 
     def test_frame(self):
@@ -57,5 +57,5 @@ def test_frame(self):
                     self.system.data["coords"][0][ii][jj],
                     ovito_posis[ii][jj],
                     places=6,
-                    msg="posis[%d][%d] failed" % (ii, jj),
+                    msg="posis[%d][%d] failed" % (ii, jj),  # noqa: UP031
                 )
diff --git a/tests/pwmat/config_ref_oh.py b/tests/pwmat/config_ref_oh.py
index ad546019a..cde0e3ec0 100644
--- a/tests/pwmat/config_ref_oh.py
+++ b/tests/pwmat/config_ref_oh.py
@@ -32,7 +32,7 @@ def test_cell(self):
                     self.system.data["cells"][0][ii][jj],
                     ovito_cell[ii][jj],
                     places=6,
-                    msg="cell[%d][%d] failed" % (ii, jj),
+                    msg="cell[%d][%d] failed" % (ii, jj),  # noqa: UP031
                 )
 
     def test_frame(self):
@@ -43,5 +43,5 @@ def test_frame(self):
                     self.system.data["coords"][0][ii][jj],
                     ovito_posis[ii][jj],
                     places=6,
-                    msg="posis[%d][%d] failed" % (ii, jj),
+                    msg="posis[%d][%d] failed" % (ii, jj),  # noqa: UP031
                 )
diff --git a/tests/pymatgen_data/deepmd/set.000/box.npy b/tests/pymatgen_data/deepmd/set.000/box.npy
new file mode 100644
index 000000000..566749c13
Binary files /dev/null and b/tests/pymatgen_data/deepmd/set.000/box.npy differ
diff --git a/tests/pymatgen_data/deepmd/set.000/coord.npy b/tests/pymatgen_data/deepmd/set.000/coord.npy
new file mode 100644
index 000000000..f30e67882
Binary files /dev/null and b/tests/pymatgen_data/deepmd/set.000/coord.npy differ
diff --git a/tests/pymatgen_data/deepmd/set.000/energy.npy b/tests/pymatgen_data/deepmd/set.000/energy.npy
new file mode 100644
index 000000000..3e591429a
Binary files /dev/null and b/tests/pymatgen_data/deepmd/set.000/energy.npy differ
diff --git a/tests/pymatgen_data/deepmd/set.000/force.npy b/tests/pymatgen_data/deepmd/set.000/force.npy
new file mode 100644
index 000000000..4546a8143
Binary files /dev/null and b/tests/pymatgen_data/deepmd/set.000/force.npy differ
diff --git a/tests/pymatgen_data/deepmd/set.000/virial.npy b/tests/pymatgen_data/deepmd/set.000/virial.npy
new file mode 100644
index 000000000..74292ab5f
Binary files /dev/null and b/tests/pymatgen_data/deepmd/set.000/virial.npy differ
diff --git a/tests/pymatgen_data/deepmd/type.raw b/tests/pymatgen_data/deepmd/type.raw
new file mode 100644
index 000000000..479b28ae4
--- /dev/null
+++ b/tests/pymatgen_data/deepmd/type.raw
@@ -0,0 +1,98 @@
+1
+1
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+3
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
diff --git a/tests/pymatgen_data/deepmd/type_map.raw b/tests/pymatgen_data/deepmd/type_map.raw
new file mode 100644
index 000000000..e2400b8fe
--- /dev/null
+++ b/tests/pymatgen_data/deepmd/type_map.raw
@@ -0,0 +1,4 @@
+Fe
+Li
+O
+P
diff --git a/tests/test_abacus_pw_scf.py b/tests/test_abacus_pw_scf.py
index 20751f819..0d2bdef5b 100644
--- a/tests/test_abacus_pw_scf.py
+++ b/tests/test_abacus_pw_scf.py
@@ -158,7 +158,7 @@ def test_noforcestress_job(self):
         # check below will not throw error
         system_ch4 = dpdata.LabeledSystem("abacus.scf", fmt="abacus/scf")
         # check the returned force is empty
-        self.assertFalse(system_ch4.data["forces"])
+        self.assertFalse(system_ch4.data["forces"].size)
         self.assertTrue("virials" not in system_ch4.data)
         # test append self
         system_ch4.append(system_ch4)
diff --git a/tests/test_abacus_spin.py b/tests/test_abacus_spin.py
index df8ef7d0d..4fcbc904c 100644
--- a/tests/test_abacus_spin.py
+++ b/tests/test_abacus_spin.py
@@ -16,6 +16,8 @@ def setUp(self):
     def tearDown(self):
         if os.path.isdir(self.dump_path):
             shutil.rmtree(self.dump_path)
+        if os.path.isfile("abacus.spin/INPUT"):
+            os.remove("abacus.spin/INPUT")
 
     def test_scf(self):
         os.system("cp abacus.spin/INPUT.scf abacus.spin/INPUT")
@@ -31,7 +33,7 @@ def test_scf(self):
             decimal=8,
         )
         np.testing.assert_almost_equal(
-            data["mag_forces"][0],
+            data["force_mags"][0],
             [
                 [-0.0000175013, -0.0000418680, -0.3669618965],
                 [-0.0000161517, -0.0000195198, -0.3669821632],
@@ -47,7 +49,40 @@ def test_scf(self):
         sys2 = dpdata.LabeledSystem(self.dump_path, fmt="deepmd/npy")
         np.testing.assert_almost_equal(data["spins"], sys2.data["spins"], decimal=8)
         np.testing.assert_almost_equal(
-            data["mag_forces"], sys2.data["mag_forces"], decimal=8
+            data["force_mags"], sys2.data["force_mags"], decimal=8
+        )
+
+    def test_scf_nspin2(self):
+        os.system("cp abacus.spin/INPUT.scf.nspin2 abacus.spin/INPUT")
+        mysys = dpdata.LabeledSystem("abacus.spin", fmt="abacus/scf")
+        data = mysys.data
+        self.assertAlmostEqual(data["energies"][0], -6818.719409466637)
+        np.testing.assert_almost_equal(
+            data["spins"][0],
+            [
+                [0, 0, 2.4000001004],
+                [0, 0, 2.3999994597],
+            ],
+            decimal=8,
+        )
+        np.testing.assert_almost_equal(
+            data["force_mags"][0],
+            [
+                [0, 0, -0.3669618965],
+                [0, 0, -0.3669821632],
+            ],
+            decimal=8,
+        )
+
+        # dump to deepmd-npy
+        mysys.to(file_name=self.dump_path, fmt="deepmd/npy")
+        self.assertTrue(os.path.isfile(f"{self.dump_path}/set.000/spin.npy"))
+        self.assertTrue(os.path.isfile(f"{self.dump_path}/set.000/force_mag.npy"))
+
+        sys2 = dpdata.LabeledSystem(self.dump_path, fmt="deepmd/npy")
+        np.testing.assert_almost_equal(data["spins"], sys2.data["spins"], decimal=8)
+        np.testing.assert_almost_equal(
+            data["force_mags"], sys2.data["force_mags"], decimal=8
         )
 
     def test_relax(self):
@@ -84,9 +119,9 @@ def test_relax(self):
             ]
         )
         self.assertEqual(len(data["spins"]), 3)
-        self.assertEqual(len(data["mag_forces"]), 3)
+        self.assertEqual(len(data["force_mags"]), 3)
         np.testing.assert_almost_equal(data["spins"], spins_ref, decimal=8)
-        np.testing.assert_almost_equal(data["mag_forces"], magforces_ref, decimal=8)
+        np.testing.assert_almost_equal(data["force_mags"], magforces_ref, decimal=8)
 
         # dump to deepmd-npy
         mysys.to(file_name=self.dump_path, fmt="deepmd/npy")
@@ -96,7 +131,7 @@ def test_relax(self):
         sys2 = dpdata.LabeledSystem(self.dump_path, fmt="deepmd/npy")
         np.testing.assert_almost_equal(data["spins"], sys2.data["spins"], decimal=8)
         np.testing.assert_almost_equal(
-            data["mag_forces"], sys2.data["mag_forces"], decimal=8
+            data["force_mags"], sys2.data["force_mags"], decimal=8
         )
 
     def test_md(self):
@@ -141,9 +176,9 @@ def test_md(self):
             ]
         )
         self.assertEqual(len(data["spins"]), 4)
-        self.assertEqual(len(data["mag_forces"]), 4)
+        self.assertEqual(len(data["force_mags"]), 4)
         np.testing.assert_almost_equal(data["spins"], spins_ref, decimal=8)
-        np.testing.assert_almost_equal(data["mag_forces"], magforces_ref, decimal=8)
+        np.testing.assert_almost_equal(data["force_mags"], magforces_ref, decimal=8)
 
         # dump to deepmd-npy
         mysys.to(file_name=self.dump_path, fmt="deepmd/npy")
@@ -153,5 +188,21 @@ def test_md(self):
         sys2 = dpdata.LabeledSystem(self.dump_path, fmt="deepmd/npy")
         np.testing.assert_almost_equal(data["spins"], sys2.data["spins"], decimal=8)
         np.testing.assert_almost_equal(
-            data["mag_forces"], sys2.data["mag_forces"], decimal=8
+            data["force_mags"], sys2.data["force_mags"], decimal=8
         )
+
+    def test_read_stru_spin(self):
+        mysys = dpdata.System("abacus.spin/STRU.spin", fmt="abacus/stru")
+        self.assertTrue("spins" in mysys.data)
+        print(mysys.data["spins"])
+
+        """
+    0.0000000000     0.000000000     0.000000000 mag 0 0 2
+    0.1000000000     0.1000000000     0.1000000000 mag 3
+    0.2000000000     0.2000000000     0.2000000000 mag 3 angle1 90
+    0.3000000000     0.3000000000     0.3000000000 mag 3 4 0 angle1 90 angle2 90
+        """
+        np.testing.assert_almost_equal(mysys.data["spins"][0][0], [0, 0, 2], decimal=8)
+        np.testing.assert_almost_equal(mysys.data["spins"][0][1], [0, 0, 3], decimal=8)
+        np.testing.assert_almost_equal(mysys.data["spins"][0][2], [3, 0, 0], decimal=8)
+        np.testing.assert_almost_equal(mysys.data["spins"][0][3], [0, 5, 0], decimal=8)
diff --git a/tests/test_abacus_stru_dump.py b/tests/test_abacus_stru_dump.py
index 084a5473c..ee0a3082e 100644
--- a/tests/test_abacus_stru_dump.py
+++ b/tests/test_abacus_stru_dump.py
@@ -29,6 +29,10 @@ def test_dump_stru(self):
         )
         myfilecmp(self, "abacus.scf/stru_test", "STRU_tmp")
 
+    def test_dump_stru_without_pporb(self):
+        self.system_ch4.to("stru", "STRU_tmp", mass=[12, 1])
+        self.assertTrue(os.path.isfile("STRU_tmp"))
+
     def test_dumpStruLinkFile(self):
         os.makedirs("abacus.scf/tmp", exist_ok=True)
         self.system_ch4.to(
@@ -51,6 +55,49 @@ def test_dumpStruLinkFile(self):
         if os.path.isdir("abacus.scf/tmp"):
             shutil.rmtree("abacus.scf/tmp")
 
+    def test_dump_stru_pporb_mismatch(self):
+        with self.assertRaises(KeyError, msg="pp_file is a dict and lack of pp for H"):
+            self.system_ch4.to(
+                "stru",
+                "STRU_tmp",
+                mass=[12, 1],
+                pp_file={"C": "C.upf", "O": "O.upf"},
+                numerical_orbital={"C": "C.orb", "H": "H.orb"},
+            )
+
+        with self.assertRaises(
+            ValueError, msg="pp_file is a list and lack of pp for H"
+        ):
+            self.system_ch4.to(
+                "stru",
+                "STRU_tmp",
+                mass=[12, 1],
+                pp_file=["C.upf"],
+                numerical_orbital={"C": "C.orb", "H": "H.orb"},
+            )
+
+        with self.assertRaises(
+            KeyError, msg="numerical_orbital is a dict and lack of orbital for H"
+        ):
+            self.system_ch4.to(
+                "stru",
+                "STRU_tmp",
+                mass=[12, 1],
+                pp_file={"C": "C.upf", "H": "H.upf"},
+                numerical_orbital={"C": "C.orb", "O": "O.orb"},
+            )
+
+        with self.assertRaises(
+            ValueError, msg="numerical_orbital is a list and lack of orbital for H"
+        ):
+            self.system_ch4.to(
+                "stru",
+                "STRU_tmp",
+                mass=[12, 1],
+                pp_file=["C.upf", "H.upf"],
+                numerical_orbital=["C.orb"],
+            )
+
     def test_dump_spinconstrain(self):
         self.system_ch4.to(
             "stru",
@@ -98,8 +145,52 @@ def test_dump_spin(self):
 4
 5.416431453540 4.011298860305 3.511161492417 1 1 1 mag 4.000000000000 5.000000000000 6.000000000000
 4.131588222365 4.706745191323 4.431136645083 1 1 1 mag 1.000000000000 1.000000000000 1.000000000000
-5.630930319126 5.521640894956 4.450356541303 1 1 1 mag 2.000000000000 2.000000000000 2.000000000000
-5.499851012568 4.003388899277 5.342621842622 1 1 1 mag 3.000000000000 3.000000000000 3.000000000000
+5.630930319126 5.521640894956 4.450356541303 0 0 0 mag 2.000000000000 2.000000000000 2.000000000000
+5.499851012568 4.003388899277 5.342621842622 0 0 0 mag 3.000000000000 3.000000000000 3.000000000000
+"""
+        self.assertTrue(stru_ref in c)
+
+    def test_dump_move_from_vasp(self):
+        self.system = dpdata.System()
+        self.system.from_vasp_poscar(os.path.join("poscars", "POSCAR.oh.c"))
+        self.system.to(
+            "abacus/stru",
+            "STRU_tmp",
+            pp_file={"O": "O.upf", "H": "H.upf"},
+        )
+        assert os.path.isfile("STRU_tmp")
+        with open("STRU_tmp") as f:
+            c = f.read()
+
+        stru_ref = """O
+0.0
+1
+0.000000000000 0.000000000000 0.000000000000 1 1 0
+H
+0.0
+1
+1.262185604418 0.701802783513 0.551388341420 0 0 0
+"""
+        self.assertTrue(stru_ref in c)
+
+        self.system.to(
+            "abacus/stru",
+            "STRU_tmp",
+            pp_file={"O": "O.upf", "H": "H.upf"},
+            move=[[True, False, True], [False, True, False]],
+        )
+        assert os.path.isfile("STRU_tmp")
+        with open("STRU_tmp") as f:
+            c = f.read()
+
+        stru_ref = """O
+0.0
+1
+0.000000000000 0.000000000000 0.000000000000 1 0 1
+H
+0.0
+1
+1.262185604418 0.701802783513 0.551388341420 0 1 0
 """
         self.assertTrue(stru_ref in c)
 
diff --git a/tests/test_custom_data_type.py b/tests/test_custom_data_type.py
index 230475ab3..08dc9eba2 100644
--- a/tests/test_custom_data_type.py
+++ b/tests/test_custom_data_type.py
@@ -120,4 +120,6 @@ def test_to_deepmd_hdf5(self):
     def test_from_deepmd_hdf5(self):
         self.system.to_deepmd_hdf5("data_bar.h5")
         x = dpdata.LabeledSystem("data_bar.h5", fmt="deepmd/hdf5")
+        print(self.system.data.keys())
+        print(x.data.keys())
         np.testing.assert_allclose(x.data["bar"], self.bar)
diff --git a/tests/test_deepmd_spin.py b/tests/test_deepmd_spin.py
new file mode 100644
index 000000000..3b5c99a38
--- /dev/null
+++ b/tests/test_deepmd_spin.py
@@ -0,0 +1,48 @@
+from __future__ import annotations
+
+import os
+import shutil
+import unittest
+
+from context import dpdata
+
+
+class TestDeepmdReadSpinNPY(unittest.TestCase):
+    def setUp(self):
+        self.tmp_save_path = "tmp.deepmd.spin/dump-tmp"
+
+    def tearDown(self):
+        if os.path.exists(self.tmp_save_path):
+            shutil.rmtree(self.tmp_save_path)
+
+    def check_Fe16(self, system):
+        self.assertTrue("spins" in system.data)
+        self.assertTrue("force_mags" in system.data)
+        self.assertEqual(system.data["spins"].shape, (2, 16, 3))
+        self.assertEqual(system.data["force_mags"].shape, (2, 16, 3))
+
+    def test_read_spin_npy(self):
+        system = dpdata.LabeledSystem("tmp.deepmd.spin/Fe16-npy", fmt="deepmd/npy")
+        self.check_Fe16(system)
+
+        system.to("deepmd/npy", self.tmp_save_path)
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.tmp_save_path, "set.000/spin.npy"))
+        )
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.tmp_save_path, "set.000/force_mag.npy"))
+        )
+
+    def test_read_spin_raw(self):
+        system = dpdata.LabeledSystem("tmp.deepmd.spin/Fe16-raw", fmt="deepmd/raw")
+        self.check_Fe16(system)
+
+        system.to("deepmd/raw", self.tmp_save_path)
+        self.assertTrue(os.path.isfile(os.path.join(self.tmp_save_path, "spin.raw")))
+        self.assertTrue(
+            os.path.isfile(os.path.join(self.tmp_save_path, "force_mag.raw"))
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_from_pymatgen.py b/tests/test_from_pymatgen.py
deleted file mode 100644
index 7689a9d5e..000000000
--- a/tests/test_from_pymatgen.py
+++ /dev/null
@@ -1,32 +0,0 @@
-from __future__ import annotations
-
-import os
-import unittest
-
-from comp_sys import CompSys
-from context import dpdata
-
-try:
-    from pymatgen.core import Structure  # noqa: F401
-
-    exist_module = True
-except Exception:
-    exist_module = False
-
-
-@unittest.skipIf(not exist_module, "skip pymatgen")
-class TestFormPytmatgen(unittest.TestCase, CompSys):
-    def setUp(self):
-        structure = Structure.from_file(os.path.join("poscars", "POSCAR.P42nmc"))
-        self.system_1 = dpdata.System(structure, fmt="pymatgen/structure")
-        self.system_2 = dpdata.System(
-            os.path.join("poscars", "POSCAR.P42nmc"), fmt="poscar"
-        )
-        self.places = 6
-        self.e_places = 6
-        self.f_places = 6
-        self.v_places = 6
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/tests/test_lammps_spin.py b/tests/test_lammps_spin.py
new file mode 100644
index 000000000..d3d58920e
--- /dev/null
+++ b/tests/test_lammps_spin.py
@@ -0,0 +1,153 @@
+from __future__ import annotations
+
+import os
+import shutil
+import unittest
+
+import numpy as np
+from context import dpdata
+
+from dpdata.lammps.dump import get_spin
+
+TRAJ_NO_ID = """ITEM: TIMESTEP
+0
+ITEM: NUMBER OF ATOMS
+17
+ITEM: BOX BOUNDS xy xz yz pp pp pp
+-4.0080511965879438e-02 5.7039029418994556e+00 -5.9179115295410201e-03
+1.4436085788922526e-02 5.6674744441011660e+00 -1.1487414836883500e-02
+7.8239288740356017e-03 5.6734038274259646e+00 6.8277359008788905e-04
+ITEM: ATOMS type x y z c_spin[1] c_spin[2] c_spin[3] c_spin[4] c_spin[5] c_spin[6] c_spin[7] c_spin[8] c_spin[9] c_spin[10]
+1 0.00141160 5.64868599 0.01005602 1.54706291 0.00000000 0.00000000 1.00000000 -1.40772100 -2.03739417 -1522.64797384 -0.00397809 -0.00190426 -0.00743976
+1 5.65283939 5.57449025 2.84281508 1.54412869 0.00000000 0.00000000 1.00000000 7.75304092 6.48949619 -1512.84926162 -0.00637234 -0.00733168 0.00661107
+1 0.00066480 2.78022036 0.01010716 1.54612979 0.00000000 0.00000000 1.00000000 -0.93618575 1.92206111 -1520.80305011 -0.00316673 0.00177893 -0.00744575
+1 5.65233666 2.85374747 2.84289453 1.54439093 0.00000000 0.00000000 1.00000000 8.11012818 -6.49922039 -1514.31557088 -0.00569217 0.00741000 0.00640353
+1 2.82063515 5.64869321 0.01007552 1.54714250 0.00000000 0.00000000 1.00000000 2.49070852 -2.14456666 -1523.53038650 0.00478410 -0.00213962 -0.00751154
+1 2.89579803 5.57439179 2.84287630 1.54415032 0.00000000 0.00000000 1.00000000 -8.03062338 6.63950296 -1513.41291897 0.00440396 -0.00717185 0.00633657
+1 2.82151287 2.78010538 0.01016303 1.54619615 0.00000000 0.00000000 1.00000000 2.71859584 1.98482729 -1521.34149633 0.00533453 0.00194532 -0.00745901
+1 2.89637049 2.85377083 2.84297332 1.54440023 0.00000000 0.00000000 1.00000000 -7.76758760 -6.67134514 -1514.43304618 0.00505040 0.00743195 0.00630302
+1 1.41106492 1.38817482 1.72302072 1.18134529 0.00000000 0.00000000 1.00000000 0.27170165 -0.00426695 -444.22843899 0.00100237 -0.00002725 -0.03385965
+ 1 1.41105247 1.38807861 3.96314606 1.18153407 0.00000000 0.00000000 1.00000000 -0.07722674 0.01368756 -337.08703133 -0.00066982 0.00007487 0.07887183
+ 1 1.41105864 4.21395432 1.43987180 1.71989299 0.00000000 0.00000000 1.00000000 -0.01511106 0.00320081 -1653.34500916 0.00010421 0.00007248 0.00634401
+ 1 1.41104843 4.21387554 4.24576823 1.71989825 0.00000000 0.00000000 1.00000000 -0.71645898 0.05923960 -1640.68070568 -0.00117959 0.00006676 -0.01467806
+ 1 4.27433865 1.38779084 1.43977211 1.72010048 0.00000000 0.00000000 1.00000000 0.45899480 0.03956420 -1653.36356942 0.00051885 0.00002313 0.00911600
+ 1 4.27436799 1.38772964 4.24586490 1.72010133 0.00000000 0.00000000 1.00000000 0.38385331 0.07301994 -1642.06086017 -0.00002034 0.00010335 -0.01688908
+ 1 4.27435427 4.21452597 1.39359689 1.65590121 0.00000000 0.00000000 1.00000000 -0.01658773 -0.06159007 -1659.12744163 0.00006470 -0.00006420 -0.01342935
+ 1 4.27434583 4.21455469 4.29208004 1.65592002 0.00000000 0.00000000 1.00000000 -0.15590720 -0.03252166 -1654.84697132 -0.00066755 -0.00003915 -0.00482188
+ 2 1.41096761 1.38958048 0.01029027 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00048351 -0.00022876 -0.00645195"""
+
+
+class TestLmp(unittest.TestCase):
+    def setUp(self):
+        self.tmp_system = dpdata.System(
+            os.path.join("poscars", "conf.lmp"), type_map=["O", "H"]
+        )
+        self.tmp_system.data["spins"] = [[[3, 4, 0], [0, 4, 3]]]
+        self.lmp_coord_name = "tmp.lmp"
+
+    def tearDown(self):
+        pass  # if os.path.isfile(self.lmp_coord_name):os.remove(self.lmp_coord_name)
+
+    def test_dump_input(self):
+        self.tmp_system.to("lammps/lmp", self.lmp_coord_name)
+        self.assertTrue(os.path.isfile(self.lmp_coord_name))
+        with open(self.lmp_coord_name) as f:
+            c = f.read()
+
+        coord_ref = """     1      1    0.0000000000    0.0000000000    0.0000000000    0.6000000000    0.8000000000    0.0000000000    5.0000000000
+     2      2    1.2621856000    0.7018028000    0.5513885000    0.0000000000    0.8000000000    0.6000000000    5.0000000000"""
+        self.assertTrue(coord_ref in c)
+
+    def test_dump_input_zero_spin(self):
+        self.tmp_system.data["spins"] = [[[0, 0, 0], [0, 0, 0]]]
+        self.tmp_system.to("lammps/lmp", self.lmp_coord_name)
+        self.assertTrue(os.path.isfile(self.lmp_coord_name))
+        with open(self.lmp_coord_name) as f:
+            c = f.read()
+        coord_ref = """     1      1    0.0000000000    0.0000000000    0.0000000000    0.0000000000    0.0000000000    1.0000000000    0.0000000000
+     2      2    1.2621856000    0.7018028000    0.5513885000    0.0000000000    0.0000000000    1.0000000000    0.0000000000"""
+        self.assertTrue(coord_ref in c)
+
+    def test_read_input(self):
+        # check if dpdata can read the spins
+        tmp_system = dpdata.System(
+            "lammps/spin.lmp", fmt="lammps/lmp", type_map=["O", "H"]
+        )
+        self.assertTrue((tmp_system.data["spins"][0] == [[3, 4, 0], [0, 4, 3]]).all())
+
+        tmp_system.to(file_name="lammps/dump", fmt="deepmd/npy")
+        self.assertTrue(os.path.isfile("lammps/dump/set.000/spin.npy"))
+
+        if os.path.isdir("lammps/dump"):
+            shutil.rmtree("lammps/dump")
+
+
+class TestDump(unittest.TestCase):
+    def test_read_dump_spin(self):
+        tmp_system = dpdata.System(
+            "lammps/traj.dump",
+            fmt="lammps/dump",
+            type_map=["O", "H"],
+            input_file="lammps/in.lmp",
+        )
+        self.assertTrue(len(tmp_system.data["spins"]) == 2)
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][0][0], [0, 0, 1.54706291], decimal=8
+        )
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][0][1], [0, 0, 1.54412869], decimal=8
+        )
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][0][-2], [0, 0, 1.65592002], decimal=8
+        )
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][0][-1], [0, 0, 0], decimal=8
+        )
+
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][1][0],
+            [0.21021514724299958, 1.0123821159859323, -0.6159960941686954],
+            decimal=8,
+        )
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][1][1],
+            [1.0057302798645609, 0.568273899191638, -0.2363447073875224],
+            decimal=8,
+        )
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][1][-2],
+            [-0.28075943761984146, -1.2845200151690905, -0.0201237855118935],
+            decimal=8,
+        )
+        np.testing.assert_almost_equal(
+            tmp_system.data["spins"][1][-1], [0, 0, 0], decimal=8
+        )
+
+        tmp_system.to(file_name="lammps/dump", fmt="deepmd/npy")
+        self.assertTrue(os.path.isfile("lammps/dump/set.000/spin.npy"))
+
+        if os.path.isdir("lammps/dump"):
+            shutil.rmtree("lammps/dump")
+
+    def test_read_dump_partial_spin(self):
+        # test if dpdata can read the spins when the spin data is not complete
+        with self.assertWarns(UserWarning) as cm:
+            tmp_system = dpdata.System(
+                "lammps/traj_partial_spin.dump",
+                fmt="lammps/dump",
+                type_map=["O", "H"],
+                input_file="lammps/in.lmp",
+            )
+            self.assertTrue("spins" not in tmp_system.data)
+
+        self.assertIn("Warning: spin info is not found in frame", str(cm.warning))
+
+    def test_get_spin_failed(self):
+        with self.assertWarns(UserWarning) as cm:
+            spin = get_spin(
+                TRAJ_NO_ID.split("\n"),
+                ["c_spin[1]", "c_spin[2]", "c_spin[3]", "c_spin[4]"],
+            )
+            self.assertTrue(spin is None)
+
+        self.assertIn("Error processing spin data:", str(cm.warning))
diff --git a/tests/test_msd.py b/tests/test_msd.py
index 7148b0b5a..5d26db645 100644
--- a/tests/test_msd.py
+++ b/tests/test_msd.py
@@ -28,6 +28,6 @@ def test_msd(self):
         # print(msd)
         ncomp = msd.shape[0]
         for ii in range(ncomp):
-            self.assertAlmostEqual(msd0[ii], ii * ii, msg="msd0[%d]" % ii)
-            self.assertAlmostEqual(msd1[ii], ii * ii * 4, msg="msd1[%d]" % ii)
-            self.assertAlmostEqual(msd[ii], (msd0[ii] + msd1[ii]) * 0.5, "msd[%d]" % ii)
+            self.assertAlmostEqual(msd0[ii], ii * ii, msg="msd0[%d]" % ii)  # noqa: UP031
+            self.assertAlmostEqual(msd1[ii], ii * ii * 4, msg="msd1[%d]" % ii)  # noqa: UP031
+            self.assertAlmostEqual(msd[ii], (msd0[ii] + msd1[ii]) * 0.5, "msd[%d]" % ii)  # noqa: UP031
diff --git a/tests/test_pymatgen_structure.py b/tests/test_pymatgen_structure.py
new file mode 100644
index 000000000..1e93829e1
--- /dev/null
+++ b/tests/test_pymatgen_structure.py
@@ -0,0 +1,61 @@
+from __future__ import annotations
+
+import os
+import unittest
+
+from comp_sys import CompSys, IsNoPBC, IsPBC
+from context import dpdata
+
+try:
+    from pymatgen.core import Structure  # noqa: F401
+
+    exist_module = True
+except Exception:
+    exist_module = False
+
+
+@unittest.skipIf(not exist_module, "skip pymatgen")
+class TestFormPytmatgen(unittest.TestCase, CompSys):
+    def setUp(self):
+        structure = Structure.from_file(os.path.join("poscars", "POSCAR.P42nmc"))
+        self.system_1 = dpdata.System(structure, fmt="pymatgen/structure")
+        self.system_2 = dpdata.System(
+            os.path.join("poscars", "POSCAR.P42nmc"), fmt="poscar"
+        )
+        self.places = 6
+        self.e_places = 6
+        self.f_places = 6
+        self.v_places = 6
+
+
+@unittest.skipIf(not exist_module, "skip pymatgen")
+class TestFormToPytmatgen(unittest.TestCase, CompSys, IsPBC):
+    def setUp(self):
+        self.system = dpdata.System("pymatgen_data/deepmd/", fmt="deepmd/npy")
+        self.system_1 = self.system
+        self.system_2 = dpdata.System().from_pymatgen_structure(
+            self.system.to("pymatgen/structure")[0]
+        )
+        self.places = 6
+        self.e_places = 6
+        self.f_places = 6
+        self.v_places = 6
+
+
+@unittest.skipIf(not exist_module, "skip pymatgen")
+class TestFormToPytmatgenNopbc(unittest.TestCase, CompSys, IsNoPBC):
+    def setUp(self):
+        self.system = dpdata.System("pymatgen_data/deepmd/", fmt="deepmd/npy")
+        self.system.data["nopbc"] = True
+        self.system_1 = self.system
+        self.system_2 = dpdata.System().from_pymatgen_structure(
+            self.system.to("pymatgen/structure")[0]
+        )
+        self.places = 6
+        self.e_places = 6
+        self.f_places = 6
+        self.v_places = 6
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_remove_pbc.py b/tests/test_remove_pbc.py
index d70a2f028..7770b6fd8 100644
--- a/tests/test_remove_pbc.py
+++ b/tests/test_remove_pbc.py
@@ -39,7 +39,7 @@ def test_remove(self):
                     self.assertAlmostEqual(
                         sys["cells"][ff][ii][jj],
                         ref[ii][jj],
-                        msg="%d %d %d" % (ff, ii, jj),
+                        msg="%d %d %d" % (ff, ii, jj),  # noqa: UP031
                     )
             dists = []
             for ii in range(sys.get_natoms()):
diff --git a/tests/test_system_apply_pbc.py b/tests/test_system_apply_pbc.py
index 2114cf6a8..22dab3baa 100644
--- a/tests/test_system_apply_pbc.py
+++ b/tests/test_system_apply_pbc.py
@@ -27,7 +27,7 @@ def test_pbc(self):
                     self.assertAlmostEqual(
                         sys["coords"][ii][jj][dd],
                         bk_coord[ii][jj][dd],
-                        msg="coord[%d][%d][%d] failed" % (ii, jj, dd),
+                        msg="coord[%d][%d][%d] failed" % (ii, jj, dd),  # noqa: UP031
                     )
 
 
diff --git a/tests/test_vasp_poscar_dump.py b/tests/test_vasp_poscar_dump.py
index 62f215986..d55228117 100644
--- a/tests/test_vasp_poscar_dump.py
+++ b/tests/test_vasp_poscar_dump.py
@@ -34,6 +34,21 @@ def setUp(self):
         self.system = dpdata.System()
         self.system.from_vasp_poscar("tmp.POSCAR")
 
+    def test_dump_move_flags(self):
+        tmp_system = dpdata.System()
+        tmp_system.from_vasp_poscar(os.path.join("poscars", "POSCAR.oh.c"))
+        tmp_system.to_vasp_poscar("tmp.POSCAR")
+        self.system = dpdata.System()
+        self.system.from_vasp_poscar("tmp.POSCAR")
+        with open("tmp.POSCAR") as f:
+            content = f.read()
+
+        stru_ref = """Cartesian
+   0.0000000000    0.0000000000    0.0000000000 T T F
+   1.2621856044    0.7018027835    0.5513883414 F F F
+"""
+        self.assertTrue(stru_ref in content)
+
 
 class TestPOSCARSkipZeroAtomNumb(unittest.TestCase):
     def tearDown(self):
diff --git a/tests/test_vasp_poscar_to_system.py b/tests/test_vasp_poscar_to_system.py
index 7457d33d2..9e5d7f37f 100644
--- a/tests/test_vasp_poscar_to_system.py
+++ b/tests/test_vasp_poscar_to_system.py
@@ -14,6 +14,37 @@ def setUp(self):
         self.system = dpdata.System()
         self.system.from_vasp_poscar(os.path.join("poscars", "POSCAR.oh.c"))
 
+    def test_move_flags(self):
+        expected = np.array([[[True, True, False], [False, False, False]]])
+        self.assertTrue(np.array_equal(self.system["move"], expected))
+
+
+class TestPOSCARMoveFlags(unittest.TestCase):
+    def setUp(self):
+        self.tmp_file = "POSCAR.tmp.1"
+
+    def tearDown(self):
+        if os.path.exists(self.tmp_file):
+            os.remove(self.tmp_file)
+
+    def test_move_flags_error1(self):
+        with self.assertRaisesRegex(RuntimeError, "Invalid move flags.*?"):
+            dpdata.System().from_vasp_poscar(os.path.join("poscars", "POSCAR.oh.err1"))
+
+    def test_move_flags_error2(self):
+        with self.assertRaisesRegex(RuntimeError, "Invalid move flag: a"):
+            dpdata.System().from_vasp_poscar(os.path.join("poscars", "POSCAR.oh.err2"))
+
+    def test_move_flags_error3(self):
+        system = dpdata.System().from_vasp_poscar(
+            os.path.join("poscars", "POSCAR.oh.c")
+        )
+        system.data["move"] = np.array([[[True, True], [False, False]]])
+        with self.assertRaisesRegex(
+            RuntimeError, "Invalid move flags:.*?should be a list of 3 bools"
+        ):
+            system.to_vasp_poscar(self.tmp_file)
+
 
 class TestPOSCARDirect(unittest.TestCase, TestPOSCARoh):
     def setUp(self):
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/box.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/box.npy
new file mode 100644
index 000000000..882bf1ff8
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/box.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/coord.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/coord.npy
new file mode 100644
index 000000000..329503e22
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/coord.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/energy.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/energy.npy
new file mode 100644
index 000000000..4bcbcf236
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/energy.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/force.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/force.npy
new file mode 100644
index 000000000..292fc6805
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/force.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/force_mag.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/force_mag.npy
new file mode 100644
index 000000000..ff0a009e8
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/force_mag.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/spin.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/spin.npy
new file mode 100644
index 000000000..ae00e0e5c
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/spin.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/set.000/virial.npy b/tests/tmp.deepmd.spin/Fe16-npy/set.000/virial.npy
new file mode 100644
index 000000000..6c2980463
Binary files /dev/null and b/tests/tmp.deepmd.spin/Fe16-npy/set.000/virial.npy differ
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/type.raw b/tests/tmp.deepmd.spin/Fe16-npy/type.raw
new file mode 100644
index 000000000..74bc717cf
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-npy/type.raw
@@ -0,0 +1,16 @@
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
diff --git a/tests/tmp.deepmd.spin/Fe16-npy/type_map.raw b/tests/tmp.deepmd.spin/Fe16-npy/type_map.raw
new file mode 100644
index 000000000..08d5832e2
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-npy/type_map.raw
@@ -0,0 +1 @@
+Fe
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/box.raw b/tests/tmp.deepmd.spin/Fe16-raw/box.raw
new file mode 100644
index 000000000..1a0f0c2b9
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/box.raw
@@ -0,0 +1,2 @@
+4.852406756199998838e+00 0.000000000000000000e+00 0.000000000000000000e+00 2.308035870000000114e-02 4.868205980599999982e+00 0.000000000000000000e+00 2.678336190000000019e-02 -1.273453610000000076e-02 7.575130651200000287e+00
+4.852406756199998838e+00 0.000000000000000000e+00 0.000000000000000000e+00 2.308035870000000114e-02 4.868205980599999982e+00 0.000000000000000000e+00 2.678336190000000019e-02 -1.273453610000000076e-02 7.575130651200000287e+00
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/coord.raw b/tests/tmp.deepmd.spin/Fe16-raw/coord.raw
new file mode 100644
index 000000000..bf290591a
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/coord.raw
@@ -0,0 +1,2 @@
+7.265197249999999496e-02 4.801689380600000057e+00 7.557116530000000054e+00 2.502718440599999816e+00 4.848670014499999681e+00 7.528592681200000136e+00 6.829729569999999639e-02 2.503290553899999882e+00 7.542453851200000337e+00 2.387861248100000111e+00 2.392921626100000143e+00 1.227113879999999942e-02 1.248625939999999934e+00 1.235871750000000047e+00 1.962114260000000110e+00 3.579928339999999931e+00 1.308363060000000022e+00 1.944950500000000027e+00 1.261952210000000019e+00 3.670674079999999950e+00 1.937314930000000102e+00 3.648521740000000069e+00 3.517678329999999853e+00 1.842547099999999993e+00 7.917075379999999329e-02 4.015369999999999877e-03 3.783381850000000046e+00 2.364675399999999872e+00 3.114458000000000157e-02 3.819137730000000008e+00 4.877458919999999587e+00 2.472410029999999814e+00 3.767453989999999919e+00 2.441658079999999842e+00 2.449273470000000064e+00 3.787900969999999923e+00 1.207819179999999992e+00 1.185759330000000000e+00 5.700741540000000107e+00 3.672023149999999792e+00 1.266887440000000087e+00 5.690759779999999601e+00 1.225844649999999980e+00 3.688144059999999946e+00 5.677029669999999584e+00 3.729146230000000006e+00 3.577686089999999819e+00 5.653438930000000084e+00
+7.262864250000000688e-02 4.800096520600000360e+00 7.556544780000000294e+00 2.501984930599999934e+00 4.849790654499999576e+00 7.529256091199999723e+00 6.993106570000000299e-02 2.504862153899999999e+00 7.542524761200000150e+00 2.386465468099999931e+00 2.393284846099999807e+00 1.243044880000000046e-02 1.250970219999999911e+00 1.235900269999999912e+00 1.961877150000000070e+00 3.578310469999999910e+00 1.309157640000000011e+00 1.944905080000000064e+00 1.263246410000000042e+00 3.669684970000000046e+00 1.937208040000000020e+00 3.646978709999999957e+00 3.518230940000000029e+00 1.843518699999999955e+00 8.087652380000000252e-02 2.412619999999999869e-03 3.782949939999999955e+00 2.363124420000000114e+00 3.159122000000000319e-02 3.818357489999999910e+00 2.538490380000000060e-02 2.473290130000000087e+00 3.767864050000000020e+00 2.441273620000000033e+00 2.447853410000000007e+00 3.785945680000000202e+00 1.206648970000000043e+00 1.183960899999999983e+00 5.700463349999999707e+00 3.673867459999999863e+00 1.268293499999999963e+00 5.690915719999999567e+00 1.223096860000000063e+00 3.690037439999999780e+00 5.677984109999999696e+00 3.731158119999999911e+00 3.576031940000000020e+00 5.654460040000000021e+00
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/energy.raw b/tests/tmp.deepmd.spin/Fe16-raw/energy.raw
new file mode 100644
index 000000000..5e49119ed
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/energy.raw
@@ -0,0 +1,2 @@
+-5.454623401600000216e+04
+-5.454620955500000127e+04
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/force.raw b/tests/tmp.deepmd.spin/Fe16-raw/force.raw
new file mode 100644
index 000000000..3a50d1ef4
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/force.raw
@@ -0,0 +1,2 @@
+-1.106855285000000050e-01 1.000727386199999902e+00 2.883449123999999864e-01 -1.051338108999999993e-01 -2.162821185999999896e-01 1.511491688999999916e-01 -6.572733814000000230e-01 -1.100767982000000034e+00 1.898751854000000072e-01 7.363416830000000246e-01 2.437835339999999962e-01 -6.857924790000000104e-02 -6.825158524000000115e-01 -7.689053999999999334e-02 -4.883741223999999992e-01 9.032410482999999601e-01 -1.615540451700000046e+00 -3.581084193999999887e-01 -4.005809030999999965e-01 5.476169250000000044e-02 -2.907007447000000133e-01 3.671111684000000208e-01 1.718198793700000060e+00 3.327191424000000008e-01 -1.116415134899999995e+00 2.404410170999999874e-01 1.111270488000000028e-01 1.152163319400000052e+00 -2.274326597999999922e-01 -1.102415283000000024e-01 -2.244907613999999885e-01 -3.542942425999999756e-01 2.837310079999999859e-02 2.335184547999999971e-01 2.195733411000000135e-01 3.262157500999999793e-01 4.219752305000000203e-01 5.337641029999999898e-01 -1.105175815000000034e-01 -4.199524229000000108e-01 -8.862139848000000208e-01 -2.222332230000000075e-02 6.769434261000000230e-01 -5.677320351000000320e-01 -1.064142056000000058e-01 -7.742465350999999663e-01 1.033904147000000107e+00 1.273548630999999931e-01
+-1.255811501999999913e-01 1.024865058300000076e+00 3.048681815000000017e-01 -9.745471190000000450e-02 -2.231689940999999899e-01 1.499777952999999941e-01 -6.818135133999999598e-01 -1.119233886199999972e+00 2.084442448000000114e-01 7.586097389999999496e-01 2.483870280000000097e-01 -7.030565599999999415e-02 -7.146105834000000545e-01 -8.766770979999999325e-02 -4.958102687000000008e-01 9.329233330999999740e-01 -1.625196258799999960e+00 -3.685665191999999979e-01 -4.263045145000000091e-01 6.338842760000000520e-02 -2.973771901000000217e-01 3.982496473000000092e-01 1.715933820300000079e+00 3.187726661999999989e-01 -1.134394543700000035e+00 2.659650671000000077e-01 1.065910659000000038e-01 1.171039510700000008e+00 -2.463114885000000087e-01 -9.709863879999999781e-02 -2.318063747999999924e-01 -3.792864573999999989e-01 2.052522660000000165e-02 2.441585086999999987e-01 2.434691978000000090e-01 3.556763707000000219e-01 4.498362350999999815e-01 5.612586234999999979e-01 -1.060362118999999964e-01 -4.527059999999999973e-01 -8.982924847000000046e-01 -2.215639380000000019e-02 7.197915827999999960e-01 -5.931029316000000495e-01 -1.250889162999999915e-01 -8.099371649000000062e-01 1.048992988399999993e+00 1.175842438999999962e-01
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/force_mag.raw b/tests/tmp.deepmd.spin/Fe16-raw/force_mag.raw
new file mode 100644
index 000000000..bc600e4e6
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/force_mag.raw
@@ -0,0 +1,2 @@
+-8.534410200000000155e-03 -1.003950009999999941e-02 -1.198012942999999997e-01 1.857438689999999989e-02 7.876794299999999999e-03 -1.158623106999999985e-01 5.790393100000000050e-03 2.902189410000000080e-02 -1.471329148999999947e-01 -1.627021049999999996e-02 -1.243863620000000038e-02 -1.431188442000000116e-01 1.983349000000000136e-03 -1.628727090000000113e-02 1.338718134000000115e-01 -1.675990100000000066e-02 2.088854290000000113e-02 1.289023514000000037e-01 7.555144599999999953e-03 2.097859640000000053e-02 1.455209802000000108e-01 -1.285195899999999917e-03 -5.159665590000000329e-02 1.251078151999999999e-01 1.242615849999999939e-02 5.394605999999999490e-04 -1.211362175000000041e-01 -1.565749859999999993e-02 4.241893000000000052e-03 -1.189902988999999994e-01 -8.991798299999999300e-03 8.944375299999999737e-03 -1.409907068000000041e-01 9.810309000000000130e-04 1.193021600000000043e-03 -1.438089848999999920e-01 -5.407705300000000166e-03 -1.619149939999999996e-02 1.402259667999999904e-01 2.398664899999999878e-03 1.248675119999999999e-02 1.249879150000000050e-01 -1.821928399999999922e-03 1.193867860000000052e-02 1.449831044000000069e-01 1.506308529999999916e-02 -2.018797589999999850e-02 1.371671825999999927e-01
+-8.948930500000000510e-03 -1.040653769999999936e-02 -1.205046651000000052e-01 1.876058129999999893e-02 8.277887200000000334e-03 -1.160853092000000003e-01 6.415928999999999903e-03 2.929689569999999882e-02 -1.481234593999999904e-01 -1.695458160000000089e-02 -1.269999410000000051e-02 -1.433629164000000122e-01 3.804538599999999875e-03 -1.687982010000000158e-02 1.344380966000000011e-01 -1.629795430000000045e-02 2.116517249999999919e-02 1.286579196000000036e-01 8.172215400000000171e-03 2.095696290000000073e-02 1.456517300999999875e-01 -1.081445799999999945e-03 -5.074005670000000107e-02 1.255218261999999918e-01 1.252434070000000080e-02 8.360950000000000239e-04 -1.213422064999999939e-01 -1.578363149999999909e-02 4.818495100000000364e-03 -1.191078504000000066e-01 -8.358033399999999224e-03 8.972258800000000589e-03 -1.413513788999999909e-01 1.125360099999999980e-03 1.306719399999999942e-03 -1.443008279999999921e-01 -5.158918000000000184e-03 -1.691780309999999996e-02 1.409441345999999928e-01 2.871428599999999859e-03 1.337089040000000000e-02 1.251760003999999871e-01 -1.211985000000000059e-03 1.288838999999999957e-02 1.463661454000000084e-01 1.788364250000000141e-02 -2.094214189999999906e-02 1.372846507000000105e-01
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/spin.raw b/tests/tmp.deepmd.spin/Fe16-raw/spin.raw
new file mode 100644
index 000000000..d231af2ae
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/spin.raw
@@ -0,0 +1,2 @@
+5.563920879999999930e-02 8.026123020000000552e-02 2.258392821700000219e+00 -2.267161628999999901e-01 -8.533636129999999653e-02 2.278170776100000072e+00 -6.958726640000000019e-02 -2.977513280999999901e-01 2.152763507299999990e+00 1.632430896000000020e-01 1.433680660999999967e-01 2.177797364699999960e+00 2.457801589999999980e-02 2.197890640000000062e-01 -2.175622943500000073e+00 2.034624560000000137e-01 -1.423908618999999987e-01 -2.227208061699999853e+00 -6.516286079999999981e-02 -2.002624396000000084e-01 -2.182219526199999926e+00 2.907600839999999881e-02 6.091792472000000513e-01 -2.112683871999999852e+00 -1.055395694000000051e-01 -4.771315269999999720e-02 2.195893971999999916e+00 1.742744213000000097e-01 -8.640444440000000226e-02 2.171807686100000190e+00 1.224203283000000059e-01 -1.428901210999999893e-01 2.129426829199999816e+00 2.808756769999999972e-02 -5.308771360000000272e-02 2.136427688399999969e+00 3.904636110000000299e-02 1.628157172999999958e-01 -2.162790827599999854e+00 -4.121785079999999712e-02 -1.434962392999999947e-01 -2.246669910100000056e+00 -1.821095639999999877e-02 -1.583821810000000108e-01 -2.097885662099999937e+00 -1.871916124000000126e-01 2.405559167000000109e-01 -2.183222719799999822e+00
+6.349788520000000658e-02 8.477241620000000588e-02 2.256907346999999842e+00 -2.246460507999999967e-01 -8.900857460000000110e-02 2.278921277699999859e+00 -7.359221209999999846e-02 -2.977365738000000150e-01 2.152393609499999805e+00 1.737881296000000075e-01 1.472407542999999985e-01 2.176163003700000154e+00 1.668498000000000030e-03 2.266339256000000080e-01 -2.172927181499999971e+00 1.942038004999999956e-01 -1.459203109999999970e-01 -2.226649295900000158e+00 -7.325918500000000444e-02 -2.010864508000000106e-01 -2.177766917800000090e+00 2.292739540000000159e-02 5.983890829000000355e-01 -2.116732264199999936e+00 -1.073844010000000043e-01 -5.292913860000000165e-02 2.194751144900000117e+00 1.736157278000000082e-01 -9.463951419999999670e-02 2.170929507399999903e+00 1.114323799999999975e-01 -1.447396732000000019e-01 2.129971659400000217e+00 2.287651989999999955e-02 -5.641734639999999856e-02 2.136625202099999843e+00 3.100547039999999865e-02 1.706691644999999979e-01 -2.161995885000000062e+00 -5.254033629999999705e-02 -1.529427707999999886e-01 -2.248103780000000107e+00 -2.916460049999999860e-02 -1.674864367999999992e-01 -2.095804000599999828e+00 -2.216720739999999967e-01 2.463171910999999992e-01 -2.178967231200000132e+00
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/type.raw b/tests/tmp.deepmd.spin/Fe16-raw/type.raw
new file mode 100644
index 000000000..74bc717cf
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/type.raw
@@ -0,0 +1,16 @@
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/type_map.raw b/tests/tmp.deepmd.spin/Fe16-raw/type_map.raw
new file mode 100644
index 000000000..08d5832e2
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/type_map.raw
@@ -0,0 +1 @@
+Fe
diff --git a/tests/tmp.deepmd.spin/Fe16-raw/virial.raw b/tests/tmp.deepmd.spin/Fe16-raw/virial.raw
new file mode 100644
index 000000000..82f12f31c
--- /dev/null
+++ b/tests/tmp.deepmd.spin/Fe16-raw/virial.raw
@@ -0,0 +1,2 @@
+7.557628675103746474e+00 -3.674192989340970272e-01 -7.797686352528614417e-01 -3.674192989340970272e-01 7.717710405205966850e+00 3.814578682575011093e-01 -7.797686352528614417e-01 3.814578682575011093e-01 7.217849924461329891e+00
+7.708199193940652449e+00 -3.671245997549927309e-01 -7.796781473872188561e-01 -3.671245997549927309e-01 7.854600132463434470e+00 3.768168653463239703e-01 -7.796781473872188561e-01 3.768168653463239703e-01 7.297379300880924902e+00
diff --git a/tests/xyz/xyz_unittest.field.xyz b/tests/xyz/xyz_unittest.field.xyz
index 9900be6b2..a0e799e04 100644
--- a/tests/xyz/xyz_unittest.field.xyz
+++ b/tests/xyz/xyz_unittest.field.xyz
@@ -16,7 +16,7 @@ energy=0.2397023e+01    virial="159.582203018    1.23282341824e-05    0.18783516
  
 
   10   
-virial=" -3.68535041825    1.63204257089e-06    -4.28008468355    1.63204257089e-06    18.1630123797    0.0    -4.28008468355    0.0    3.03073243091 "    Lattice=" 9.217080809   0.0   0.0   4.86e-07   6.431267224   0.0   4.201562981   4.1e-08   2.205334915 "  Properties='species:S:1:pos:R:3:Z:I:1:force:R:3 ' energy='-58.342497      '
+virial=" -3.68535041825    1.63204257089e-06    -4.28008468355    1.63204257089e-06    18.1630123797    0.0    -4.28008468355    0.0    3.03073243091 "    Lattice=" 9.217080809   0.0   0.0   4.86e-07   6.431267224   0.0   4.201562981   4.1e-08   2.205334915 "  Properties='species:S:1:pos:R:3:Z:I:1:force:R:3 ' energy='-5.834249877929687500e+01      '
 B   12.38023000000   3.21563000000   2.13103000000   5    -9.35551000000  -0.00000000000  -0.22364600000
 C   12.96783000000   4.32879000000   2.14172000000   6     7.05653600000  11.19171000000   1.80087100000
 C   12.96783000000   2.10247000000   2.14172000000   6     7.05653500000 -11.19171000000   1.80087100000
@@ -28,7 +28,7 @@ C    2.71900000000   1.71197000000   0.17210000000   6    -3.43041200000   0.056
 C    2.71900000000   4.71930000000   0.17210000000   6    -3.43041200000  -0.05665700000  -0.44035700000
 C    9.69981000000   3.21563000000   0.55395000000   6    -2.18778000000  -0.00000000000  -1.00760300000
 10
-energy=-56.397425    virial="-5.61261501333    -6.34809444383e-07    6.34809444383e-07    -6.34809444383e-07    -5.61261501333    6.34809444383e-07    6.34809444383e-07    6.34809444383e-07    4.91130356636"    Lattice="5.512073672   0.0   0.0   4.16e-07   5.512073672   0.0   2.53e-07   2.53e-07   3.347159197"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
+energy=-5.639742660522460938e+01    virial="-5.61261501333    -6.34809444383e-07    6.34809444383e-07    -6.34809444383e-07    -5.61261501333    6.34809444383e-07    6.34809444383e-07    6.34809444383e-07    4.91130356636"    Lattice="5.512073672   0.0   0.0   4.16e-07   5.512073672   0.0   2.53e-07   2.53e-07   3.347159197"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
 B    4.08581000000   4.57623000000   2.94238000000   5     0.03989300000  -1.87257900000   0.97202400000
 C    1.42627000000   0.93584000000   2.94238000000   5    -0.03989300000   1.87257900000   0.97202400000
 C    0.93584000000   4.08581000000   2.94238000000   5     1.87257900000   0.03989300000   0.97202400000
diff --git a/tests/xyz/xyz_unittest.sort.xyz b/tests/xyz/xyz_unittest.sort.xyz
index cbd15c274..17b14eb60 100644
--- a/tests/xyz/xyz_unittest.sort.xyz
+++ b/tests/xyz/xyz_unittest.sort.xyz
@@ -13,7 +13,7 @@ B    3.19567000000   7.85111000000   1.52250000000   5    -7.25565700000   8.829
 C    2.27564000000   0.00000000000   1.03309000000   6    -0.00000000000  -0.00000000000   0.00000000000
 B    2.27564000000   4.41538000000   1.03309000000   5    -0.00000000000  -0.00000000000   0.00000000000
 10
-energy=-58.342497    virial="-3.68535041825    1.63204257089e-06    -4.28008468355    1.63204257089e-06    18.1630123797    0.0    -4.28008468355    0.0    3.03073243091"    Lattice="9.217080809   0.0   0.0   4.86e-07   6.431267224   0.0   4.201562981   4.1e-08   2.205334915"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
+energy=-5.834249877929687500e+01    virial="-3.68535041825    1.63204257089e-06    -4.28008468355    1.63204257089e-06    18.1630123797    0.0    -4.28008468355    0.0    3.03073243091"    Lattice="9.217080809   0.0   0.0   4.86e-07   6.431267224   0.0   4.201562981   4.1e-08   2.205334915"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
 C   12.96783000000   4.32879000000   2.14172000000   6     7.05653600000  11.19171000000   1.80087100000
 C   12.96783000000   2.10247000000   2.14172000000   6     7.05653500000 -11.19171000000   1.80087100000
 C    7.95424000000   1.03482000000   2.10290000000   6     1.71010500000  -3.80357800000  -0.12402100000
@@ -25,7 +25,7 @@ C    2.71900000000   1.71197000000   0.17210000000   6    -3.43041200000   0.056
 C    2.71900000000   4.71930000000   0.17210000000   6    -3.43041200000  -0.05665700000  -0.44035700000
 C    9.69981000000   3.21563000000   0.55395000000   6    -2.18778000000  -0.00000000000  -1.00760300000
 10
-energy=-56.397425    virial="-5.61261501333    -6.34809444383e-07    6.34809444383e-07    -6.34809444383e-07    -5.61261501333    6.34809444383e-07    6.34809444383e-07    6.34809444383e-07    4.91130356636"    Lattice="5.512073672   0.0   0.0   4.16e-07   5.512073672   0.0   2.53e-07   2.53e-07   3.347159197"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
+energy=-5.639742660522460938e+01    virial="-5.61261501333    -6.34809444383e-07    6.34809444383e-07    -6.34809444383e-07    -5.61261501333    6.34809444383e-07    6.34809444383e-07    6.34809444383e-07    4.91130356636"    Lattice="5.512073672   0.0   0.0   4.16e-07   5.512073672   0.0   2.53e-07   2.53e-07   3.347159197"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
 B    4.08581000000   4.57623000000   2.94238000000   5     0.03989300000  -1.87257900000   0.97202400000
 C    1.42627000000   0.93584000000   2.94238000000   5    -0.03989300000   1.87257900000   0.97202400000
 C    0.93584000000   4.08581000000   2.94238000000   5     1.87257900000   0.03989300000   0.97202400000
diff --git a/tests/xyz/xyz_unittest.xyz b/tests/xyz/xyz_unittest.xyz
index e179e3e1f..28b6ab31c 100644
--- a/tests/xyz/xyz_unittest.xyz
+++ b/tests/xyz/xyz_unittest.xyz
@@ -13,7 +13,7 @@ C    2.27564000000   1.56002000000   1.03309000000   6    -0.00000000000  -4.519
 C    2.27564000000   7.27074000000   1.03309000000   6     0.00000000000   4.51935800000   0.00000000000
 C    2.27564000000   0.00000000000   1.03309000000   6    -0.00000000000  -0.00000000000   0.00000000000
 10
-energy=-58.342497    virial="-3.68535041825    1.63204257089e-06    -4.28008468355    1.63204257089e-06    18.1630123797    0.0    -4.28008468355    0.0    3.03073243091"    Lattice="9.217080809   0.0   0.0   4.86e-07   6.431267224   0.0   4.201562981   4.1e-08   2.205334915"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
+energy=-5.834249877929687500e+01    virial="-3.68535041825    1.63204257089e-06    -4.28008468355    1.63204257089e-06    18.1630123797    0.0    -4.28008468355    0.0    3.03073243091"    Lattice="9.217080809   0.0   0.0   4.86e-07   6.431267224   0.0   4.201562981   4.1e-08   2.205334915"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
 B   12.38023000000   3.21563000000   2.13103000000   5    -9.35551000000  -0.00000000000  -0.22364600000
 C   12.96783000000   4.32879000000   2.14172000000   6     7.05653600000  11.19171000000   1.80087100000
 C   12.96783000000   2.10247000000   2.14172000000   6     7.05653500000 -11.19171000000   1.80087100000
@@ -25,7 +25,7 @@ C    2.71900000000   1.71197000000   0.17210000000   6    -3.43041200000   0.056
 C    2.71900000000   4.71930000000   0.17210000000   6    -3.43041200000  -0.05665700000  -0.44035700000
 C    9.69981000000   3.21563000000   0.55395000000   6    -2.18778000000  -0.00000000000  -1.00760300000
 10
-energy=-56.397425    virial="-5.61261501333    -6.34809444383e-07    6.34809444383e-07    -6.34809444383e-07    -5.61261501333    6.34809444383e-07    6.34809444383e-07    6.34809444383e-07    4.91130356636"    Lattice="5.512073672   0.0   0.0   4.16e-07   5.512073672   0.0   2.53e-07   2.53e-07   3.347159197"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
+energy=-5.639742660522460938e+01    virial="-5.61261501333    -6.34809444383e-07    6.34809444383e-07    -6.34809444383e-07    -5.61261501333    6.34809444383e-07    6.34809444383e-07    6.34809444383e-07    4.91130356636"    Lattice="5.512073672   0.0   0.0   4.16e-07   5.512073672   0.0   2.53e-07   2.53e-07   3.347159197"     Properties=species:S:1:pos:R:3:Z:I:1:force:R:3
 B    4.08581000000   4.57623000000   2.94238000000   5     0.03989300000  -1.87257900000   0.97202400000
 C    1.42627000000   0.93584000000   2.94238000000   5    -0.03989300000   1.87257900000   0.97202400000
 C    0.93584000000   4.08581000000   2.94238000000   5     1.87257900000   0.03989300000   0.97202400000