diff --git a/dpdata/system.py b/dpdata/system.py
index 9ba6ad840..b588bf898 100644
--- a/dpdata/system.py
+++ b/dpdata/system.py
@@ -434,6 +434,35 @@ def apply_pbc(self) :
         ncoord = ncoord % 1
         self.data['coords'] = np.matmul(ncoord, self.data['cells'])
 
+
+    def remove_pbc(self, protect_layer = 9):
+        """
+        This method does NOT delete the definition of the cells, it
+        (1) revises the cell to a cubic cell and ensures that the cell 
+        boundary to any atom in the system is no less than `protect_layer`
+        (2) translates the system such that the center-of-geometry of the system 
+        locates at the center of the cell.
+
+        Parameters
+        ----------
+        protect_layer : the protect layer between the atoms and the cell
+                        boundary
+        """
+        nframes = self.get_nframes()
+        natoms = self.get_natoms()
+        assert(protect_layer >= 0), "the protect_layer should be no less than 0"
+        for ff in range(nframes):
+            tmpcoord = self.data['coords'][ff]
+            cog = np.average(tmpcoord, axis = 0)
+            dist = tmpcoord - np.tile(cog, [natoms, 1])
+            max_dist = np.max(np.linalg.norm(dist, axis = 1))
+            h_cell_size = max_dist + protect_layer            
+            cell_size = h_cell_size * 2
+            shift = np.array([1,1,1]) * h_cell_size - cog
+            self.data['coords'][ff] = self.data['coords'][ff] + np.tile(shift, [natoms, 1])
+            self.data['cells'][ff] = cell_size * np.eye(3)            
+        
+
     @register_from_funcs.register_funcs("lmp")
     @register_from_funcs.register_funcs("lammps/lmp")
     def from_lammps_lmp (self, file_name, type_map = None) :
diff --git a/tests/test_remove_pbc.py b/tests/test_remove_pbc.py
new file mode 100644
index 000000000..ea2fc2202
--- /dev/null
+++ b/tests/test_remove_pbc.py
@@ -0,0 +1,37 @@
+import os
+import numpy as np
+import unittest
+from context import dpdata
+
+class TestRemovePBC(unittest.TestCase):
+
+    def test_remove(self):
+        coords = np.array([[[-1, -1, 2], [-1,-1,-3], [-1,-1, 7]], 
+                           [[ 3, -1, 3], [-1,-1, 3], [ 7,-1, 3]]], dtype = float)
+        cogs = np.average(coords, axis = 1)        
+        data = {'atom_names' : ['A', 'B'],
+                'atom_numbs' : [1, 2],
+                'atom_types' : np.array([1, 0, 1], dtype = int),
+                'orig': np.array([0, 0, 0]),
+                'coords': coords,
+                'cells': np.random.random([2, 3, 3]),
+        }
+        sys = dpdata.System(data = data)
+        proct = 9.0
+        
+        mol_size = np.array([5, 4], dtype = float)
+        cell_size = (mol_size + proct) * 2.0
+
+        sys.remove_pbc(proct)
+
+        for ff in range(2):
+            ref = cell_size[ff] * np.eye(3)
+            for ii in range(3):
+                for jj in range(3):
+                    self.assertAlmostEqual(sys['cells'][ff][ii][jj], ref[ii][jj], msg = '%d %d %d' %(ff, ii, jj))
+            dists = []
+            for ii in range(sys.get_natoms()):
+                for jj in range(3):
+                    dists.append(np.abs(sys['coords'][ff][ii][jj]))
+                    dists.append(np.abs(sys['cells'][ff][jj][jj] - sys['coords'][ff][ii][jj]))
+            self.assertAlmostEqual(np.min(dists), proct)