diff --git a/README.md b/README.md
index ae041803..63da9a0c 100644
--- a/README.md
+++ b/README.md
@@ -131,26 +131,26 @@ Available properties are (nframe: number of frames in the system, natoms: total
 ## Dump data
 The data stored in `System` or `LabeledSystem` can be dumped in 'lammps/lmp' or 'vasp/poscar' format, for example:
 ```python
-d_outcar.to_lammps_lmp('conf.lmp', frame_idx=0)
+d_outcar.to('lammps/lmp', 'conf.lmp', frame_idx=0)
 ```
 The first frames of `d_outcar` will be dumped to 'conf.lmp'
 ```python
-d_outcar.to_vasp_poscar('POSCAR', frame_idx=-1)
+d_outcar.to('vasp/poscar', 'POSCAR', frame_idx=-1)
 ```
 The last frames of `d_outcar` will be dumped to 'POSCAR'.
 
 
 The data stored in `LabeledSystem` can be dumped to deepmd-kit raw format, for example
 ```python
-d_outcar.to_deepmd_raw('dpmd_raw')
+d_outcar.to('deepmd/raw', 'dpmd_raw')
 ```
 Or a simpler command:
 ```python
-dpdata.LabeledSystem('OUTCAR').to_deepmd_raw('dpmd_raw')
+dpdata.LabeledSystem('OUTCAR').to('deepmd/raw', 'dpmd_raw')
 ```
 Frame selection can be implemented by
 ```python
-dpdata.LabeledSystem('OUTCAR').sub_system([0,-1]).to_deepmd_raw('dpmd_raw')
+dpdata.LabeledSystem('OUTCAR').sub_system([0,-1]).to('deepmd/raw', 'dpmd_raw')
 ```
 by which only the first and last frames are dumped to `dpmd_raw`.
 
diff --git a/dpdata/qe/traj.py b/dpdata/qe/traj.py
index d88ea872..e9800330 100644
--- a/dpdata/qe/traj.py
+++ b/dpdata/qe/traj.py
@@ -1,7 +1,7 @@
 #!/usr/bin/python3 
 
 import numpy as np
-import dpdata
+import dpdata,warnings
 
 ry2ev = 13.605693009
 hartree2ev = 27.211386018
@@ -36,8 +36,10 @@ def convert_celldm(ibrav, celldm) :
         return celldm[0] * 0.5 * np.array([[1,1,1], [-1,1,1], [-1,-1,1]])
     elif ibrav == -3 :
         return celldm[0] * 0.5 * np.array([[-1,1,1], [1,-1,1], [1,1,-1]])
-    else :
-        raise RuntimeError('unsupported ibrav ' + str(ibrav))    
+    else :        
+        warnings.warn('unsupported ibrav ' + str(ibrav) + ' if no .cel file, the cell convertion may be wrong. ')
+        return np.eye(3)
+        #raise RuntimeError('unsupported ibrav ' + str(ibrav))
 
 def load_cell_parameters(lines) :
     blk = load_block(lines, 'CELL_PARAMETERS', 3)
@@ -98,15 +100,16 @@ def _load_pos_block(fp, natoms) :
     head = fp.readline()
     if not head:
         # print('get None')
-        return None
+        return None, None
     else :
+        ss = head.split()[0]
         blk = []
         for ii in range(natoms) :
             newline = fp.readline()
             if not newline :
-                return None
+                return None, None
             blk.append([float(jj) for jj in newline.split()])
-        return blk
+        return blk, ss
 
 
 def load_data(fname, 
@@ -115,18 +118,20 @@ def load_data(fname,
               step = 1,
               convert = 1.) :
     coords = []
+    steps = []
     cc = 0
     with open(fname) as fp:
         while True:
-            blk = _load_pos_block(fp, natoms)
+            blk, ss = _load_pos_block(fp, natoms)
             if blk == None :
                 break
             else :
                 if cc >= begin and (cc - begin) % step == 0 :
                     coords.append(blk)
+                    steps.append(ss)
             cc += 1
     coords= convert * np.array(coords)
-    return coords
+    return coords, steps
 
 
 # def load_pos(fname, natoms) :
@@ -145,14 +150,19 @@ def load_data(fname,
 
 def load_energy(fname, begin = 0, step = 1) :
     data = np.loadtxt(fname)
+    steps = []
+    for ii in data[begin::step,0]:
+        steps.append('%d'%ii)
     with open(fname) as fp:
-        line = fp.readline()
-        if line :
-            nw = len(line.split())
-        else :
-            return None
+        while True:
+            line = fp.readline()
+            if not line :
+                return None
+            if line.split()[0][0] != '#':
+                nw = len(line.split())
+                break
     data = np.reshape(data, [-1, nw])
-    return energy_convert * data[begin::step,5]
+    return energy_convert * data[begin::step,5], steps
 
 
 # def load_force(fname, natoms) :
@@ -176,35 +186,48 @@ def to_system_data(input_name, prefix, begin = 0, step = 1) :
         data['atom_types'], \
         cell \
         = load_param_file(input_name)
-    data['coords'] \
+    data['coords'], csteps\
         = load_data(prefix + '.pos', 
                     np.sum(data['atom_numbs']), 
                     begin = begin, 
                     step = step, 
                     convert = length_convert)
     data['orig'] = np.zeros(3)
-    data['cells'] = np.tile(cell, (data['coords'].shape[0], 1, 1))
-    return data
+    try :
+        data['cells'], tmp_steps \
+            = load_data(prefix + '.cel', 
+                        3, 
+                        begin = begin, 
+                        step = step, 
+                        convert = length_convert)
+        assert(csteps == tmp_steps), "the step key between files are not consistent"
+    except FileNotFoundError :
+        data['cells'] = np.tile(cell, (data['coords'].shape[0], 1, 1))
+    return data, csteps
 
 
 def to_system_label(input_name, prefix, begin = 0, step = 1) :
     atom_names, atom_numbs, atom_types, cell = load_param_file(input_name)
-    energy = load_energy(prefix + '.evp', 
-                         begin = begin, 
-                         step = step)
-    force = load_data(prefix + '.for', 
-                      np.sum(atom_numbs), 
-                      begin = begin,
-                      step = step, 
-                      convert = force_convert)
-    return energy, force
+    energy, esteps = load_energy(prefix + '.evp', 
+                                 begin = begin, 
+                                 step = step)
+    force, fsteps = load_data(prefix + '.for', 
+                              np.sum(atom_numbs), 
+                              begin = begin,
+                              step = step, 
+                              convert = force_convert)
+    assert(esteps == fsteps), "the step key between files are not consistent "
+    return energy, force, esteps
 
 
 if __name__ == '__main__':
-    atom_names, atom_numbs, atom_types, cell = load_param_file('oh-md.in')
-    coords = load_pos('oh-md.pos', np.sum(atom_numbs))
+    prefix='nacl'
+    atom_names, atom_numbs, atom_types, cell = load_param_file(prefix+'.in')
+    coords = load_data(prefix+'.pos', np.sum(atom_numbs))
+    cells = load_data(prefix+'.cel', 3)
     print(atom_names)
     print(atom_numbs)
     print(atom_types)
-    print(cell)
+    print(cells)
     print(coords.shape)
+    print(cells.shape)
diff --git a/dpdata/system.py b/dpdata/system.py
index c129b877..b588bf89 100644
--- a/dpdata/system.py
+++ b/dpdata/system.py
@@ -1,6 +1,7 @@
 #%%
 import os
 import glob
+import inspect
 import numpy as np
 import dpdata.lammps.lmp
 import dpdata.lammps.dump
@@ -23,6 +24,18 @@
 from dpdata.periodic_table import Element
 from dpdata.xyz.quip_gap_xyz import QuipGapxyzSystems
 
+
+class Register:
+    def __init__(self):
+        self.funcs = {}
+
+    def register_funcs(self, fmt):
+        def decorator(func):
+            self.funcs[fmt] = func
+            return func
+        return decorator
+
+
 class System (MSONable) :
     '''
     The data System
@@ -97,30 +110,37 @@ def __init__ (self,
             return
         if file_name is None :
             return
-        if fmt == 'auto':
-            fmt = os.path.basename(file_name).split('.')[-1]
-        if fmt == 'lmp' or fmt == 'lammps/lmp' :
-            self.from_lammps_lmp(file_name, type_map = type_map)
-        elif fmt == 'dump' or fmt == 'lammps/dump' :
-            self.from_lammps_dump(file_name, type_map = type_map, begin = begin, step = step)
-        elif fmt.lower() == 'poscar' or fmt.lower() == 'contcar' or fmt.lower() == 'vasp/poscar' or fmt.lower() == 'vasp/contcar':
-            self.from_vasp_poscar(file_name)
-        elif fmt == 'deepmd' or fmt == 'deepmd/raw':
-            self.from_deepmd_raw(file_name, type_map = type_map)
-        elif fmt == 'deepmd/npy':
-            self.from_deepmd_comp(file_name, type_map = type_map)
-        elif fmt == 'qe/cp/traj':
-            self.from_qe_cp_traj(file_name, begin = begin, step = step)
-        elif fmt.lower() == 'siesta/output':
-            self.from_siesta_output(file_name)
-        elif fmt.lower() == 'siesta/aimd_output':
-            self.from_siesta_aiMD_output(file_name)
-        else :
-            raise RuntimeError('unknow data format ' + fmt)
+        self.from_fmt(file_name, fmt, type_map=type_map, begin= begin, step=step)
 
         if type_map is not None:
             self.apply_type_map(type_map)
 
+    register_from_funcs = Register()
+    register_to_funcs = Register()
+
+    def from_fmt(self, file_name, fmt='auto', **kwargs):
+        fmt = fmt.lower()
+        if fmt == 'auto':
+            fmt = os.path.basename(file_name).split('.')[-1].lower()
+        from_funcs = self.register_from_funcs.funcs
+        if fmt in from_funcs:
+            func = from_funcs[fmt]
+            args = inspect.getfullargspec(func).args
+            kwargs = {kk: kwargs[kk] for kk in kwargs if kk in args}
+            func(self, file_name, **kwargs)
+        else :
+            raise RuntimeError('unknow data format ' + fmt)
+    
+    def to(self, fmt, *args, **kwargs):
+        fmt = fmt.lower()
+        to_funcs = self.register_to_funcs.funcs
+        if fmt in to_funcs:
+            func = to_funcs[fmt]
+            func_args = inspect.getfullargspec(func).args
+            kwargs = {kk: kwargs[kk] for kk in kwargs if kk in func_args}
+            func(self, *args, **kwargs)
+        else :
+            raise RuntimeError('unknow data format %s. Accepted format:' % (fmt, " ".join(to_funcs)))
 
     def __repr__(self):
         return self.__str__()
@@ -207,7 +227,7 @@ def map_atom_types(self,type_map=None):
 
         return new_atom_types
 
-
+    @register_to_funcs.register_funcs("list")
     def to_list(self):
         """
         convert system to list, usefull for data collection
@@ -415,12 +435,43 @@ def apply_pbc(self) :
         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) :
         with open(file_name) as fp:
             lines = [line.rstrip('\n') for line in fp]
             self.data = dpdata.lammps.lmp.to_system_data(lines, type_map)
         self._shift_orig_zero()
 
+    @register_to_funcs.register_funcs("pymatgen/structure")
     def to_pymatgen_structure(self):
         '''
         convert System to Pymatgen Structure obj
@@ -440,6 +491,7 @@ def to_pymatgen_structure(self):
             structures.append(structure)
         return structures
 
+    @register_to_funcs.register_funcs("ase/structure")
     def to_ase_structure(self):
         '''
         convert System to ASE Atom obj
@@ -459,6 +511,7 @@ def to_ase_structure(self):
             structures.append(structure)
         return structures
 
+    @register_to_funcs.register_funcs("lammps/lmp")
     def to_lammps_lmp(self, file_name, frame_idx = 0) :
         """
         Dump the system in lammps data format
@@ -475,7 +528,8 @@ def to_lammps_lmp(self, file_name, frame_idx = 0) :
         with open(file_name, 'w') as fp:
             fp.write(w_str)
 
-
+    @register_from_funcs.register_funcs('dump')
+    @register_from_funcs.register_funcs('lammps/dump')
     def from_lammps_dump (self,
                           file_name,
                           type_map = None,
@@ -485,13 +539,17 @@ def from_lammps_dump (self,
         self.data = dpdata.lammps.dump.system_data(lines, type_map)
         self._shift_orig_zero()
 
-
+    @register_from_funcs.register_funcs('poscar')
+    @register_from_funcs.register_funcs('contcar')
+    @register_from_funcs.register_funcs('vasp/poscar')
+    @register_from_funcs.register_funcs('vasp/contcar')
     def from_vasp_poscar(self, file_name) :
         with open(file_name) as fp:
             lines = [line.rstrip('\n') for line in fp]
             self.data = dpdata.vasp.poscar.to_system_data(lines)
         self.rot_lower_triangular()
 
+    @register_to_funcs.register_funcs("vasp/string")
     def to_vasp_string(self, frame_idx=0):
         """
         Dump the system in vasp POSCAR format string
@@ -505,6 +563,7 @@ def to_vasp_string(self, frame_idx=0):
         w_str = dpdata.vasp.poscar.from_system_data(self.data, frame_idx)
         return w_str
     
+    @register_to_funcs.register_funcs("vasp/poscar")
     def to_vasp_poscar(self, file_name, frame_idx = 0) :
         """
         Dump the system in vasp POSCAR format
@@ -520,27 +579,30 @@ def to_vasp_poscar(self, file_name, frame_idx = 0) :
         with open(file_name, 'w') as fp:
             fp.write(w_str)
 
-
+    @register_from_funcs.register_funcs('qe/cp/traj')
     def from_qe_cp_traj(self,
                         prefix,
                         begin = 0,
                         step = 1) :
-        self.data = dpdata.qe.traj.to_system_data(prefix + '.in', prefix, begin = begin, step = step)
+        self.data, _ = dpdata.qe.traj.to_system_data(prefix + '.in', prefix, begin = begin, step = step)
         self.data['coords'] \
             = dpdata.md.pbc.apply_pbc(self.data['coords'],
                                       self.data['cells'],
             )
         self.rot_lower_triangular()
 
-
+    @register_from_funcs.register_funcs('deepmd/npy')
     def from_deepmd_comp(self, folder, type_map = None) :
         self.data = dpdata.deepmd.comp.to_system_data(folder, type_map = type_map, labels = False)
 
+    @register_from_funcs.register_funcs('deepmd')
+    @register_from_funcs.register_funcs('deepmd/raw')
     def from_deepmd_raw(self, folder, type_map = None) :
         tmp_data = dpdata.deepmd.raw.to_system_data(folder, type_map = type_map, labels = False)
         if tmp_data is not None :
             self.data = tmp_data
 
+    @register_to_funcs.register_funcs("deepmd/npy")
     def to_deepmd_npy(self, folder, set_size = 5000, prec=np.float32) :
         """
         Dump the system in deepmd compressed format (numpy binary) to `folder`.
@@ -563,12 +625,14 @@ def to_deepmd_npy(self, folder, set_size = 5000, prec=np.float32) :
                                 set_size = set_size,
                                 comp_prec = prec)
 
+    @register_to_funcs.register_funcs("deepmd/raw")
     def to_deepmd_raw(self, folder) :
         """
         Dump the system in deepmd raw format to `folder`
         """
         dpdata.deepmd.raw.dump(folder, self.data) 
 
+    @register_from_funcs.register_funcs('siesta/output')
     def from_siesta_output(self, fname):
         self.data['atom_names'], \
         self.data['atom_numbs'], \
@@ -579,6 +643,7 @@ def from_siesta_output(self, fname):
             = dpdata.siesta.output.obtain_frame(fname)
         # self.rot_lower_triangular()
 
+    @register_from_funcs.register_funcs('siesta/aimd_output')
     def from_siesta_aiMD_output(self, fname):
         self.data['atom_names'], \
         self.data['atom_numbs'], \
@@ -737,6 +802,14 @@ def nopbc(self):
             return True
         return False
 
+    def shuffle(self):
+        """Shuffle frames randomly."""
+        idx = np.random.permutation(self.get_nframes())
+        for ii in ['cells', 'coords']:
+            self.data[ii] = self.data[ii][idx]
+        return idx
+
+
 def get_cell_perturb_matrix(cell_pert_fraction):
     if cell_pert_fraction<0:
         raise RuntimeError('cell_pert_fraction can not be negative')
@@ -836,38 +909,11 @@ def __init__ (self,
            return
         if file_name is None :
             return
-        if fmt == 'auto':
-            fmt = os.path.basename(file_name).split('.')[-1]
-        if fmt == 'xml' or fmt == 'XML' or fmt == 'vasp/xml' :
-            self.from_vasp_xml(file_name, begin = begin, step = step)
-        elif fmt == 'outcar' or fmt == 'OUTCAR' or fmt == 'vasp/outcar' :
-            self.from_vasp_outcar(file_name, begin = begin, step = step)
-        elif fmt == 'deepmd' or fmt == 'deepmd/raw':
-            self.from_deepmd_raw(file_name, type_map = type_map)
-        elif fmt == 'deepmd/npy':
-            self.from_deepmd_comp(file_name, type_map = type_map)
-        elif fmt == 'qe/cp/traj':
-            self.from_qe_cp_traj(file_name, begin = begin, step = step)
-        elif fmt == 'qe/pw/scf':
-            self.from_qe_pw_scf(file_name)
-        elif fmt.lower() == 'siesta/output':
-            self.from_siesta_output(file_name)
-        elif fmt.lower() == 'siesta/aimd_output':
-            self.from_siesta_aiMD_output(file_name)
-        elif fmt == 'gaussian/log':
-            self.from_gaussian_log(file_name)
-        elif fmt == 'gaussian/md':
-            self.from_gaussian_log(file_name, md=True)
-        elif fmt == 'cp2k/output':
-            self.from_cp2k_output(file_name)
-        elif fmt == 'cp2k/aimd_output':
-            self.from_cp2k_aimd_output(file_dir=file_name)
-        else :
-            raise RuntimeError('unknow data format ' + fmt)
-
+        self.from_fmt(file_name, fmt, type_map=type_map, begin= begin, step=step)
         if type_map is not None:
             self.apply_type_map(type_map)
 
+    register_from_funcs = Register()
 
     def __repr__(self):
         return self.__str__()
@@ -905,7 +951,7 @@ def has_virial(self) :
         # return ('virials' in self.data) and (len(self.data['virials']) > 0)
         return ('virials' in self.data)
 
-
+    @register_from_funcs.register_funcs('cp2k/aimd_output')
     def from_cp2k_aimd_output(self, file_dir):
         xyz_file=glob.glob("{}/*pos*.xyz".format(file_dir))[0]
         log_file=glob.glob("{}/*.log".format(file_dir))[0]
@@ -913,6 +959,8 @@ def from_cp2k_aimd_output(self, file_dir):
             l = LabeledSystem(data=info_dict)
             self.append(l)
 
+    @register_from_funcs.register_funcs('xml')
+    @register_from_funcs.register_funcs('vasp/xml')
     def from_vasp_xml(self, file_name, begin = 0, step = 1) :
         self.data['atom_names'], \
             self.data['atom_types'], \
@@ -937,7 +985,8 @@ def from_vasp_xml(self, file_name, begin = 0, step = 1) :
         # rotate the system to lammps convention
         self.rot_lower_triangular()
 
-
+    @register_from_funcs.register_funcs('outcar')
+    @register_from_funcs.register_funcs('vasp/outcar')
     def from_vasp_outcar(self, file_name, begin = 0, step = 1) :
         # with open(file_name) as fp:
         #     lines = [line.rstrip('\n') for line in fp]
@@ -979,27 +1028,30 @@ def rot_frame_lower_triangular(self, f_idx = 0) :
         self.affine_map_fv(trans, f_idx = f_idx)
         return trans
 
-
+    @register_from_funcs.register_funcs('deepmd/npy')
     def from_deepmd_comp(self, folder, type_map = None) :
         self.data = dpdata.deepmd.comp.to_system_data(folder, type_map = type_map, labels = True)
 
-
+    @register_from_funcs.register_funcs('deepmd')
+    @register_from_funcs.register_funcs('deepmd/raw')
     def from_deepmd_raw(self, folder, type_map = None) :
         tmp_data = dpdata.deepmd.raw.to_system_data(folder, type_map = type_map, labels = True)
         if tmp_data is not None :
             self.data = tmp_data
 
-
+    @register_from_funcs.register_funcs('qe/cp/traj')
     def from_qe_cp_traj(self, prefix, begin = 0, step = 1) :
-        self.data = dpdata.qe.traj.to_system_data(prefix + '.in', prefix, begin = begin, step = step)
+        self.data, cs = dpdata.qe.traj.to_system_data(prefix + '.in', prefix, begin = begin, step = step)
         self.data['coords'] \
             = dpdata.md.pbc.apply_pbc(self.data['coords'],
                                       self.data['cells'],
             )
-        self.data['energies'], self.data['forces'] \
+        self.data['energies'], self.data['forces'], es \
             = dpdata.qe.traj.to_system_label(prefix + '.in', prefix, begin = begin, step = step)
+        assert(cs == es), "the step key between files are not consistent" 
         self.rot_lower_triangular()
 
+    @register_from_funcs.register_funcs('qe/pw/scf')
     def from_qe_pw_scf(self, file_name) :
         self.data['atom_names'], \
             self.data['atom_numbs'], \
@@ -1012,6 +1064,7 @@ def from_qe_pw_scf(self, file_name) :
             = dpdata.qe.scf.get_frame(file_name)
         self.rot_lower_triangular()
     
+    @register_from_funcs.register_funcs('siesta/output')
     def from_siesta_output(self, file_name) :
         self.data['atom_names'], \
         self.data['atom_numbs'], \
@@ -1024,6 +1077,7 @@ def from_siesta_output(self, file_name) :
             = dpdata.siesta.output.obtain_frame(file_name)
         # self.rot_lower_triangular()
 
+    @register_from_funcs.register_funcs('siesta/aimd_output')
     def from_siesta_aiMD_output(self, file_name):
         self.data['atom_names'], \
         self.data['atom_numbs'], \
@@ -1035,14 +1089,19 @@ def from_siesta_aiMD_output(self, file_name):
         self.data['virials'] \
             = dpdata.siesta.aiMD_output.get_aiMD_frame(file_name)
     
+    @register_from_funcs.register_funcs('gaussian/log')
     def from_gaussian_log(self, file_name, md=False):
         try:
             self.data = dpdata.gaussian.log.to_system_data(file_name, md=md)
         except AssertionError:
             self.data['energies'], self.data['forces']= [], []
             self.data['nopbc'] = True
+    
+    @register_from_funcs.register_funcs('gaussian/md')
+    def from_gaussian_md(self, file_name):
+        self.from_gaussian_log(file_name, md=True)
 
-
+    @register_from_funcs.register_funcs('cp2k/output')
     def from_cp2k_output(self, file_name) :
         self.data['atom_names'], \
             self.data['atom_numbs'], \
@@ -1109,6 +1168,14 @@ def sort_atom_types(self):
             if ii in self.data:
                 self.data[ii] = self.data[ii][:, idx]
 
+    def shuffle(self):
+        """Also shuffle labeled data e.g. energies and forces."""
+        idx = System.shuffle(self)
+        for ii in ['energies', 'forces', 'virials', 'atom_pref']:
+            if ii in self.data:
+                self.data[ii] = self.data[ii][idx]
+        return idx
+
 
 class MultiSystems:
     '''A set containing several systems.'''
diff --git a/tests/qe.traj/traj6.cel b/tests/qe.traj/traj6.cel
new file mode 100644
index 00000000..bcba5c37
--- /dev/null
+++ b/tests/qe.traj/traj6.cel
@@ -0,0 +1,24 @@
+195
+5.157620138324029213e+00 0 0
+8.882263018337732685e-02 5.782217839848955876e+00 0
+4.531494312205264219e-02 1.857154025688012577e-01 5.855553673164934025e+00
+200
+5.359985500701728967e+00 0 0
+3.585941820098031974e-01 5.317218997480877896e+00 0
+7.606780476053129902e-01 7.811107228901693622e-01 5.715864930517207121e+00
+201
+5.235484852416917079e+00 0 0
+5.612267823951659906e-01 5.640178212424978632e+00 0
+1.561235654197479228e-02 2.819721738692628765e-01 5.868536463888631260e+00
+202
+5.735200354152383717e+00 0 0
+8.818112752306928037e-02 5.003946763988647461e+00 0
+5.346029912365874992e-01 4.453576669882056693e-01 5.311832768820492490e+00
+203
+5.969007839117177916e+00 0 0
+7.974158979687617776e-02 5.086270090760550922e+00 0
+2.286579398249665163e-01 1.449088120223311904e-01 5.588188387987865546e+00
+204
+5.308510801020571712e+00 0 0
+3.076052782312116429e-01 5.279388982187173340e+00 0
+4.321921336152507731e-01 8.121110815096156399e-01 5.301664983741235737e+00
diff --git a/tests/qe.traj/traj6.evp b/tests/qe.traj/traj6.evp
index 0717b0f2..b122f248 100644
--- a/tests/qe.traj/traj6.evp
+++ b/tests/qe.traj/traj6.evp
@@ -1,6 +1,7 @@
+# comments here
     195  9.433649E-03  1.290500E-05  0.000000E+00  1.057182E-02      -1100.03076389      -1100.03076389      -1100.03075434      -1100.03073209  1.300606E+04      -3.31610
     200  9.675537E-03  7.090171E-05  0.000000E+00  1.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
-    200  9.675537E-03  7.090171E-05  0.000000E+00  2.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
-    200  9.675537E-03  7.090171E-05  0.000000E+00  3.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
-    200  9.675537E-03  7.090171E-05  0.000000E+00  4.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
-    200  9.675537E-03  7.090171E-05  0.000000E+00  5.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
+    201  9.675537E-03  7.090171E-05  0.000000E+00  2.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
+    202  9.675537E-03  7.090171E-05  0.000000E+00  3.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
+    203  9.675537E-03  7.090171E-05  0.000000E+00  4.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
+    204  9.675537E-03  7.090171E-05  0.000000E+00  5.257319E+00      -1100.03189493      -1100.03189493      -1100.03076015      -1100.03068120  1.300606E+04      -3.31683
diff --git a/tests/qe.traj/traj6.for b/tests/qe.traj/traj6.for
index d9154efe..94d7110a 100644
--- a/tests/qe.traj/traj6.for
+++ b/tests/qe.traj/traj6.for
@@ -4,15 +4,15 @@
     200  0.00967554
      0.13452693616244E-01    -0.68546285592109E-02    -0.28341150176630E-02
      0.23859341498989E-01    -0.15249791118496E-01    -0.77171937738584E-02
-    200  0.00967554
+ 201  0.00967554
      2.13452693616244E-01    -0.68546285592109E-02    -0.28341150176630E-02
      0.23859341498989E-01    -0.15249791118496E-01    -0.77171937738584E-02
-    200  0.00967554
+ 202  0.00967554
      3.13452693616244E-01    -0.68546285592109E-02    -0.28341150176630E-02
      0.23859341498989E-01    -0.15249791118496E-01    -0.77171937738584E-02
-    200  0.00967554
+ 203  0.00967554
      4.13452693616244E-01    -0.68546285592109E-02    -0.28341150176630E-02
      0.23859341498989E-01    -0.15249791118496E-01    -0.77171937738584E-02
-    200  0.00967554
+ 204  0.00967554
      5.13452693616244E-01    -0.68546285592109E-02    -0.28341150176630E-02
      0.23859341498989E-01    -0.15249791118496E-01    -0.77171937738584E-02
diff --git a/tests/qe.traj/traj6.pos b/tests/qe.traj/traj6.pos
index 19af926f..965e8c08 100644
--- a/tests/qe.traj/traj6.pos
+++ b/tests/qe.traj/traj6.pos
@@ -4,15 +4,15 @@
     200  0.00967554
      0.32121855340601E+02     0.17164728967776E+02    -0.36182628967521E+01
      0.14321036136090E+01     0.54322965963498E+01     0.16038534464162E+02
-    200  0.00967554
+    201  0.00967554
      2.32121855340601E+02     0.17164728967776E+02    -0.36182628967521E+01
      0.14321036136090E+01     0.54322965963498E+01     0.16038534464162E+02
-    200  0.00967554
+    202  0.00967554
      3.32121855340601E+02     0.17164728967776E+02    -0.36182628967521E+01
      0.14321036136090E+01     0.54322965963498E+01     0.16038534464162E+02
-    200  0.00967554
+    203  0.00967554
      4.32121855340601E+02     0.17164728967776E+02    -0.36182628967521E+01
      0.14321036136090E+01     0.54322965963498E+01     0.16038534464162E+02
-    200  0.00967554
+    204  0.00967554
      5.32121855340601E+02     0.17164728967776E+02    -0.36182628967521E+01
      0.14321036136090E+01     0.54322965963498E+01     0.16038534464162E+02
diff --git a/tests/test_lammps_lmp_dump.py b/tests/test_lammps_lmp_dump.py
index d20f2984..f3354648 100644
--- a/tests/test_lammps_lmp_dump.py
+++ b/tests/test_lammps_lmp_dump.py
@@ -14,6 +14,15 @@ def setUp(self):
         self.system.from_lammps_lmp('tmp.lmp', 
                                     type_map = ['O', 'H'])
 
+class TestToFunc(unittest.TestCase, TestPOSCARoh):
+    
+    def setUp(self): 
+        tmp_system = dpdata.System(os.path.join('poscars', 'conf.lmp'), 
+                                   type_map = ['O', 'H'])
+        tmp_system.to('lammps/lmp', 'tmp.lmp')
+        self.system = dpdata.System()
+        self.system.from_fmt('tmp.lmp', fmt='lammps/lmp',
+                             type_map = ['O', 'H'])
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_qe_cp_traj.py b/tests/test_qe_cp_traj.py
index 6fab9824..b98a47a4 100644
--- a/tests/test_qe_cp_traj.py
+++ b/tests/test_qe_cp_traj.py
@@ -4,7 +4,7 @@
 from context import dpdata
 
 
-class TestPWSCFProps :
+class TestCPTRAJProps :
     def test_atom_names(self) :
         self.assertEqual(self.system.data['atom_names'], ['O','H'])
 
@@ -44,18 +44,25 @@ def test_coord(self) :
                 self.assertAlmostEqual(self.system['coords'][-1][ii][jj], coords[ii][jj])
 
 
-class TestPWSCFTraj(unittest.TestCase, TestPWSCFProps):    
+class TestCPTRAJTraj(unittest.TestCase, TestCPTRAJProps):    
 
     def setUp(self): 
         self.system = dpdata.System('qe.traj/oh-md', fmt = 'qe/cp/traj')
 
 
-class TestPWSCFLabeledTraj(unittest.TestCase, TestPWSCFProps):    
+class TestCPTRAJLabeledTraj(unittest.TestCase, TestCPTRAJProps):    
 
     def setUp(self): 
         self.system = dpdata.LabeledSystem('qe.traj/oh-md', fmt = 'qe/cp/traj')
 
 
+class TestConverCellDim(unittest.TestCase):    
+    def test_case_null(self):
+        cell = dpdata.qe.traj.convert_celldm(8, [1, 1, 1])
+        ref = np.eye(3)
+        for ii in range(3):
+            for jj in range(3):
+                self.assertAlmostEqual(cell[ii][jj], ref[ii][jj])
 
 
 if __name__ == '__main__':
diff --git a/tests/test_qe_cp_traj_skipload.py b/tests/test_qe_cp_traj_skipload.py
index dc02e632..c87baa31 100644
--- a/tests/test_qe_cp_traj_skipload.py
+++ b/tests/test_qe_cp_traj_skipload.py
@@ -38,3 +38,21 @@ def setUp(self):
         self.f_places = 6
         self.v_places = 4
 
+    def test_cell(self):
+        ref_cell = [5.359985500701728967e+00, 0, 0,
+                    3.585941820098031974e-01, 5.317218997480877896e+00, 0,
+                    7.606780476053129902e-01, 7.811107228901693622e-01, 5.715864930517207121e+00 ]
+        ref_cell = 0.52917721067 * np.array(ref_cell).reshape(3,3)
+        
+        for ii in range(3) :
+            for jj in range(3) :
+                self.assertEqual(self.system_1.data['cells'][0][ii][jj], ref_cell[ii][jj])
+
+        ref_cell = [5.308510801020571712e+00, 0, 0,
+                    3.076052782312116429e-01, 5.279388982187173340e+00, 0,
+                    4.321921336152507731e-01, 8.121110815096156399e-01, 5.301664983741235737e+00]
+        ref_cell = 0.52917721067 * np.array(ref_cell).reshape(3,3)
+        
+        for ii in range(3) :
+            for jj in range(3) :
+                self.assertEqual(self.system_1.data['cells'][-1][ii][jj], ref_cell[ii][jj])
diff --git a/tests/test_qe_pw_scf.py b/tests/test_qe_pw_scf.py
index eb1c96bb..391e763e 100644
--- a/tests/test_qe_pw_scf.py
+++ b/tests/test_qe_pw_scf.py
@@ -38,6 +38,7 @@ def test_cell(self) :
         for ii in range(cell.shape[0]) :
             for jj in range(cell.shape[1]) :
                 self.assertAlmostEqual(self.system_h2o.data['cells'][0][ii][jj], cell[ii][jj])
+        fp.close()
 
 
     def test_coord(self) :
diff --git a/tests/test_remove_pbc.py b/tests/test_remove_pbc.py
new file mode 100644
index 00000000..ea2fc220
--- /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)
diff --git a/tests/test_shuffle.py b/tests/test_shuffle.py
new file mode 100644
index 00000000..aa2d7a7e
--- /dev/null
+++ b/tests/test_shuffle.py
@@ -0,0 +1,21 @@
+import unittest
+from context import dpdata
+from comp_sys import CompLabeledSys, IsPBC
+
+class TestDeepmdLoadRaw(unittest.TestCase, CompLabeledSys, IsPBC):
+    def setUp (self) :
+        original_system = dpdata.LabeledSystem('poscars/OUTCAR.h2o.md', 
+                                             fmt = 'vasp/outcar')
+        original_system += original_system
+        original_system += original_system
+        original_system += original_system
+        self.system_1 = dpdata.LabeledSystem()
+        self.system_2 = original_system.copy()
+        idx = self.system_2.shuffle()
+        for ii in idx:
+            self.system_1.append(original_system.sub_system(ii))
+
+        self.places = 6
+        self.e_places = 6
+        self.f_places = 6
+        self.v_places = 6