Repository: deepmodeling/rid-kit (language: python)

Issue:
Dflow
1. Add labeling MD statistical analysis to remove the mean forces with high std.
2. Add CV output in MCMC projection.
3. Add training entrypoints.
4. Update docs and README.
5. Add exploration entrypoints.
6. Update MCMC plot levels.
7. Release of the dflow refactor version, 0.5.2.


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pipy_release.yml b/.github/workflows/pipy_release.yml
new file mode 100644
index 00000000..2ab7e767
--- /dev/null
+++ b/.github/workflows/pipy_release.yml
@@ -0,0 +1,34 @@
+name: Publish Python 🐍 distributions 📦 to PyPI
+
+on:
+  push:
+    tags:
+     - '*'
+
+jobs:
+  build-n-publish:
+    name: Build and publish Python 🐍 distributions 📦 to PyPI
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@master
+    - name: Set up Python 3.9
+      uses: actions/setup-python@v3
+      with:
+        python-version: '3.9'
+    - name: Install pypa/setuptools
+      run: >-
+        python -m
+        pip install wheel
+    - name: Extract tag name
+      id: tag
+      run: echo ::set-output name=TAG_NAME::$(echo $GITHUB_REF | cut -d / -f 3)
+    - name: Update version in setup.py
+      run: >-
+        sed -i "/use_scm_version=/d" setup.py && sed -i "s/{{VERSION_PLACEHOLDER}}/${{ steps.tag.outputs.TAG_NAME }}/g" setup.py
+    - name: Build a binary wheel
+      run: >-
+        python setup.py sdist bdist_wheel
+    - name: Publish distribution 📦 to PyPI
+      uses: pypa/gh-action-pypi-publish@master
+      with:
+        password: ${{ secrets.PYPI_API_TOKEN }}
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 8bb4cf74..a80efaae 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,10 +1,11 @@
 __pycache__
 *.pyc
 .nfs*
-rid.egg-info
+rid_kit.egg-info
 *.egg
 dist
 build
 tmp*
 _version.py
-.coverage
\ No newline at end of file
+.coverage
+.DS_Store
\ No newline at end of file
diff --git a/README.md b/README.md
index 3da5e824..8555e4ee 100644
--- a/README.md
+++ b/README.md
@@ -1,25 +1,32 @@
 # Table of contents
 - [About Rid-kit](#about-rid-kit)
-- [Quick Start](#quick-start)
-- [Run the workflow without k8s enviroment](#run-the-workflow-without-k8s-environment)
+- [Use Rid-kit with Bohrium k8s](#use-rid-kit-with-bohrium-k8s)
+  - [Quick Start](#quick-start)
+- [Use Rid-kit with Local k8s](#use-rid-kit-with-local-k8s)
+- [Use Rid-kit without k8s environment](#use-rid-kit-without-k8s-environment)
 - [Main procedure of RiD](#main-procedure-of-rid)
-- [Use Rid-kit](#use-rid-kit)
-- [Installation of enviroment](#installation-of-enviroment)
-- [Installation of DeepMD potential](#installation-of-dp)
+- [Preparing files in the input directory](#preparing-files-in-the-input-directory)
 - [Configure simulations](#configure-simulations)
 - [Configure machine resources](#configure-machine-resources)
+- [Download files](#download-files-from-the-workflow)
+- [Configure MCMC dimension reduction](#configure-mcmc)
+- [Installation of environment on Slurm](#installation-of-environment-on-slurm)
+- [Installation with DeepMD potential support](#installation-with-deepmd-potential-support)
 - [Troubleshooting](#troubleshooting)
 
 # About Rid-kit
 Rid-kit is a package written in Python, designed to do enhanced sampling using reinforced dynamics. It aims to learn the free energy surface on the fly during MD run, and uses it as the bias potential during the next MD run. Its advantage is the ability to use a large number of CVs (100), thus can be used to simulate conformational changes of big molecules such as in the problem of protein folding.
 
+Rid-kit is based on [dflow](https://github.com/deepmodeling/dflow), one can decide whether to use a k8s environment to run a workflow. It is recommended to use k8s environment to run the workflow unless it is not doable anyway. 
+
 For more information, check the [documentation](https://rid-kit-dflow.readthedocs.io/en/latest/).
 
-# Quick start
-Rid-kit is based on [dflow](https://github.com/deepmodeling/dflow), one can decide whether to use a k8s environment to run a workflow. It is recommended to use k8s enviroment to run the workflow unless it is not doable anyway. The dflow team provide a community version of k8s [deepmodeling k8s](https://workflows.deepmodeling.com/), making the use of Rid-kit very convenient. To use the community version of k8s, one first need to register a Bohrium account in [Bohrium](https://bohrium.dp.tech/login) and learn a few concepts (job, jobgroup, project id) in the Bohrium website documents. Then the use of rid-kit is very easy.
+# Use Rid-kit with Bohrium k8s
+## Quick start
+The dflow team provide a community version of k8s [deepmodeling k8s](https://workflows.deepmodeling.com/), making the use of Rid-kit very convenient. To use the community version of k8s, one first need to register a Bohrium account in [Bohrium](https://bohrium.dp.tech/login) and learn a few concepts (job, jobgroup, project id) in the Bohrium website documents. Then the use of rid-kit is very easy.
 
-## Set the enviroment variables
-Just set the enviroment variables based on your personal Bohrium account information by
+### Set the environment variables
+Just set the environment variables based on your personal Bohrium account information by
 
 ```
 export DFLOW_HOST=https://workflows.deepmodeling.com
@@ -31,43 +38,100 @@ export BOHRIUM_PASSWORD="<bohrium-password>"
 export BOHRIUM_PROJECT_ID="<bohrium-project-id>"
 ```
 
-## Install Rid-kit
+### Install Rid-kit
 Install the latest rid-kit
-```
-git clone git@github.com:PKUfjh/rid-kit.git
-cd rid-kit
-git checkout dflow
+
+```bash
 pip install setuptools_scm
-pip install .
+pip install -U rid-kit
 ```
-## Run an example
+
+### Run an example
 Change to the rid-kit directory
-```
+
+```bash
 cd rid-kit
 ```
-Run a example of Ala-dipeptide using dihedral as CVs (change to your own Bohrium account information)
-```
+
+Run a example of Ala-dipeptide on Bohrium using dihedral as CVs (change to your own Bohrium account information)
+
+```bash
 rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_bohrium_k8s.json
 ```
-You can also run the example on a Slurm machine (But you need to configure a conda enviroment on the slurm, see below)
-```
+
+You can also run the example on a Slurm machine (But you need to configure a conda environment on the slurm, see [Installation](docs/source/install.md))
+
+```bash
 rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_slurm_k8s.json
 ```
-You can also run the example using distance as CVs (This will use constrained md as the mean force calculator)
+
+You can specify the workflow name by providing WORKFLOW_ID after "-d", for example:
+
+```bash
+rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_bohrium_k8s.json -d ala-dipeptide-1
+```
+
+**Note that the defined workflow-id should only contain lower case alphanumeric character, and specifal character "-".**
+
+You can also specify other types of CVs such as distance or any customized CVs, for detailed explanation you can check 
+- [configure simulations](docs/source/rid_configuration.md)
+
+### Continue from an old workflow
+
+Using `resubmit` to continue from an old workflow
+
+```bash
+# suppose the original workflow id is OLD_ID
+rid resubmit -i ./tests/data/000 -c ./rid/template/rid_gmx_dis.json -m ./rid/template/machine_bohrium_k8s.json OLD_ID -d NEW_ID 
+```
+
+If you want to resubmit from a particular `iteration` and `step`:
+
+```bash
+rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json OLD_ID -t ITERATION-ID -p STEP-KEY -d NEW_ID
+```
+
+ The `ITERATION-ID`(start from 1) is just `n`th iteration the workflow has been executed. The `STEP-KEY` in rid includes the following steps: `prep-exploration`, `run-exploration`, `prep-select`, `run-select`, `prep-label`, `run-label`, `label-stats`, `collect-data`, `merge-data`, `train`, `model-devi`.
+
+### Download files from the workflow
+
+```bash
+rid download WORKFLOW_ID -p STEP-KEY -f FILE_NAME -a ITERATION_START -e ITERATION_END -o OUTPUT_DIR
 ```
-rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dis.json -m ./rid/template/machine_bohrium_k8s.json
+
+typically we want the trajectories information from each exploration step, suppose we run a workflow for 20 iterations.
+
+```bash
+rid download WORKFLOW_ID -p run-exploration -f trajectory -a 1 -e 20 -o my_protein_out
+```
+
+Detained explanation for other files in the Rid run can be found as follows:
+- [Download files](docs/source/rid_download.md)
+
+### Reduce the dimension of free energy model
+After the Rid-kit Run, the workflow will generate several numbers of free energy models (.pb). The Rid-kit currently support MCMC to reduce the dimension of the free energy model, for example:
+```bash
+rid redim -i ./test/data/models -c ./rid/template/rid_mcmc_cv_dih.json -m ./rid/template/machine_bohrium_k8s_mcmc.json
 ```
+Then you will get the projected free energy surface for ala-dipeptide
+![image1](docs/pics/mcmc_ala_1cv.png) | ![image2](docs/pics/mcmc_ala_2cv.png)
+---|---
+
+You can also include .out file representing the CV output information inside the directory specified by `-i` parameter, this will plot the CV output upon the free energy surface.
 
-Note that if you want to use constrained MD as the mean force calculator, apart from setting `method` to be `constrained` in the `label_config`, you should add `[ constraints ]` line corresponding to the `[ moleculartype ]` in your input `topology` file yourself, since gromacs specifies constraints information for each `[ moleculartype ]`.
+# Use Rid-kit with Local k8s
 
-Also note that, since gromacs only supports constrained MD for `distance CV`, the constrained MD simulation in rid-kit only supports `distance CV` at this moment.
+A tutorial on using Rid-kit with k8s environment configured by your own can be found as follows:
 
-# Run the Workflow without k8s environment
+- [Tutorial](docs/source/tutorial.md)
+
+# Use Rid-kit without k8s environment
 To run the workflow without k8s environment, one can use the `Debug` mode of `Dflow`. In this mode however, one can not monitor the workflow in the `Argo` UI.
 ## Run an example
-If one wants to run the workflow on the `Slurm` machine locally, type
-```
-DFLOW_DEBUG=1 rid submit -i /tests/data/000 -c /rid/template/rid_gmx_dih.json -m /rid/template/machine_slurm_local.json
+If one wants to run the workflow on the `Slurm` machine locally, change to the rid-kit directory and type (change to your slurm configuration)
+
+```bash
+DFLOW_DEBUG=1 rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_slurm_local.json -d ala-dipeptide-1
 ```
 
 # Main procedure of RiD
@@ -92,30 +156,30 @@ A fully connected NN will be trained via sampling data. This network will genera
 
 A more detailed description of RiD is published now, please see:
 
->  [1]  Zhang, L., Wang, H., E, W.. Reinforced dynamics for enhanced sampling in large atomic and molecular systems[J]. The Journal of chemical physics, 2018, 148(12): 124113.
+> [1]  Zhang, L., Wang, H., E, W.. Reinforced dynamics for enhanced sampling in large atomic and molecular systems[J]. The Journal of chemical physics, 2018, 148(12): 124113.
 >  
->  [2]  Wang, D., Wang, Y., Chang, J., Zhang, L., & Wang, H. (2022). Efficient sampling of high-dimensional free energy landscapes using adaptive reinforced dynamics. Nature Computational Science, 2(1), 20-29.
+> [2]  Wang, D., Wang, Y., Chang, J., Zhang, L., & Wang, H. (2022). Efficient sampling of high-dimensional free energy landscapes using adaptive reinforced dynamics. Nature Computational Science, 2(1), 20-29.
 
-# Use Rid-kit
+# Preparing files in the input directory
+Rules for preparing files in the input directory can be found in
+- [Preparing files](docs/source/prepare_files.md)
 
-A tutorial on using Rid-kit can be found as follows:
+# Configure simulations
+- [configure simulations](docs/source/rid_configuration.md)
 
-- [Tutorial](/docs/source/tutorial.ipynb)
+# Configure machine resources
+- [configure machines](docs/source/rid_machine.md)
 
-# Enviroment settings
+# Configure mcmc 
+- [configure mcmc dimension reduction](docs/source/mcmc_configuration.md)
 
-Installation of the computation enviroment can be found in
+# Installation of environment on Slurm
+If you sumit the RiD workflow to Bohrium, you do not need to install the enviroment yourself, rather Bohrium will pull the docker images to do the computation. But if you submit the workflow to Slurm machine, you will have to install the computation environment, details of the installation can be found in
 - [Installation](docs/source/install.md)
 
-# Installation-of-dp
-Installation of the DeepMD potential support can be found in
+# Installation with DeepMD potential support
+Installation of the computation environment with DeepMD potential support on slurm machines can be found in
 - [Installation of dp](docs/source/install_dp.md)
-  
-# Configure simulations
-- [configure simulations](docs/source/rid_configuration.md)
-
-# Configure machine resources
-- [configure machines](docs/source/rid_machine.md)
 
 # Workflow Synopsis
 
diff --git a/docs/pics/mcmc_ala_1cv.png b/docs/pics/mcmc_ala_1cv.png
new file mode 100644
index 00000000..4450bf18
Binary files /dev/null and b/docs/pics/mcmc_ala_1cv.png differ
diff --git a/docs/pics/mcmc_ala_2cv.png b/docs/pics/mcmc_ala_2cv.png
new file mode 100644
index 00000000..0b2f017b
Binary files /dev/null and b/docs/pics/mcmc_ala_2cv.png differ
diff --git a/docs/requirements.txt b/docs/requirements.txt
index ea1ac232..4701a2d6 100755
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -2,10 +2,13 @@ pydflow
 numpy
 scikit-learn
 mdtraj
-tensorflow
+dpdata
+matplotlib
+# tensorflow
 ipython
 jinja2==3.0.1
 docutils==0.17.1
 myst-parser
 m2r2
-nbsphinx
\ No newline at end of file
+nbsphinx
+sphinx_rtd_theme
\ No newline at end of file
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 47962629..d97cb311 100755
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -25,6 +25,8 @@
 # The full version, including alpha/beta/rc tags
 release = '2.0'
 
+nbsphinx_allow_errors = True
+
 
 # -- General configuration ---------------------------------------------------
 
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 4ba2409e..d3404dd3 100755
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -20,13 +20,16 @@ Contents
 .. toctree::
    :maxdepth: 2
 
-   Tutorial <tutorial>
-   Environment settings  <install>
-   Installation with DP potential <install_dp>
+   Use Rid-kit with Local k8s <tutorial>
+   Prepare files for simulations <prepare_files>
    Configure simulations <rid_configuration>
    Configure machine resources <rid_machine>
-   Rid API <api/api>
+   Download files from simulation <rid_download>
+   Configure mcmc dimension reduction <mcmc_configuration>
+   Installation of environment on Slurm  <install>
+   Installation with DeepMD potential support <install_dp>
    Trouble Shooting <troubleshooting/index> 
+   Rid API <api/api>
 
 Indices and tables
 ==================
diff --git a/docs/source/install.md b/docs/source/install.md
index 40c897ed..5de40126 100644
--- a/docs/source/install.md
+++ b/docs/source/install.md
@@ -51,15 +51,27 @@ export PLUMED_KERNEL=${CONDA_PREFIX}/lib/libplumedKernel.@SOEXT@
 ### Install gromacs 2022.4
 Note that in order to compile gromacs with GPU support, you need to have CUDA toolkit installed. In slurm enviroment, if the software is managed by `Enviroment Module`, one can load CUDA toolkits such as `module load cuda/11.4`. Also note that in most slurm environment, there is no network connection at computing node, but if you want to compile gromacs without fftw installed, you will have to compile it on login-in node where network connection is on.
 
+However if the network is not on in the slurm machine, you can first compile fftw manually by
+```bash
+cd fftw-3.3.8
+./configure --prefix=$CONDA_PREFIX
+make
+make install
+```
+
 ```bash
 tar -xzvf gromacs-2022.4.tar.gz
 cd gromacs-2022.4
 plumed patch -p # make sure you use the correct gromacs version
 mkdir build
 cd build
+# if the network is on
 cmake .. -DGMX_BUILD_OWN_FFTW=ON \
                 -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX \
                 -DGMX_GPU=CUDA
+# if the network is off
+cmake .. -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX \
+            -DGMX_GPU=CUDA
 make -j 4
 make install
 ```
diff --git a/docs/source/mcmc_configuration.md b/docs/source/mcmc_configuration.md
new file mode 100644
index 00000000..5a6aa234
--- /dev/null
+++ b/docs/source/mcmc_configuration.md
@@ -0,0 +1,85 @@
+# Setting Parameters of MCMC dimension reduction
+
+`RiD-kit` uses a JSON-format file to configure MCMC dimension reduction. Here we explain these parameters one by one.
+
+```JSON
+{
+        "name": "rid-mcmc",
+        "init_models": ["model_000.pb","model_001.pb","model_002.pb","model_003.pb"],
+        
+        "MCMC_Config": {
+            "cv_dimension": 2,
+            "numb_steps": 10000,
+            "numb_walkers": 500,
+            "temperature": 300,
+            "proj_info": {
+                "proj_mode": "cv",
+                "proj_cv_index": [0,1]
+            },
+            "cv_type": "dih",
+            "bins": 101
+        }
+    
+},
+{
+        "name": "rid-mcmc",
+        "init_models": ["model_000.pb","model_001.pb","model_002.pb","model_003.pb"],
+        
+        "MCMC_Config": {
+            "cv_dimension": 20,
+            "numb_steps": 100000,
+            "numb_walkers": 2000,
+            "cv_upper_bound": [4.5,5.0,6.0,4.2,4.0,4.0,3.5,4.0,4.0,3.5,3.2,3.0,3.5,3.5,4.0,4.5,4.5,3.75,3.8,4.0],
+            "cv_lower_bound": [1.0,2.0,3.5,2.8,2.5,2.0,1.2,1.5,1.5,1.0,1.8,0.5,1.0,1.0,1.5,1.5,1.5,1.75,2.2,2.6],
+            "proj_info": {
+                "proj_mode": "cv",
+                "proj_cv_index": [0,1]
+            },
+            "cv_type": "dis",
+            "bins": 101
+        }
+    
+}
+```
+* **`"init_models"`** `(List[str])` List of trained free energy models generated in the RiD run.
+* **`"cv_dimension"`** `int` An integer to represent the number of CVs used in RiD.
+* **`"numb_steps"`** `int` An integer to represent the number of steps for MCMC run.
+* **`"numb_walkers"`** `int` Number of parallel walkers used in the MCMC run.
+* **`"temperature"`**  `int` Temperature used in the RiD run and MCMC simulation.
+* **`"cv_upper_bound"`** `Optional[(List[float])]` Upper bound of the initial value in MCMC for each dimension of CVs, this is usually set for distance CVs.
+* **`"cv_lower_bound"`** `Optional[(List[float])]` Lower bound of the initial value in MCMC for each dimension of CVs, this is usually set for distance CVs.
+* **`"proj_mode"`** `str` The mode for projecting 2D CV, support `cv` and `path` mode. `cv` mode representing the common projection on the selected CV index, `path` mode representing the projection on the user defined path CV.
+* **`"proj_cv_index"`** `(List[int])` Projected index for the 2D CV. Note that 1D projection is done for all CVs.
+* **`"cv_type"`** `str` CV type for projection, support `dih` and `dis`, representing dihedral CVs and distance CVs.
+* **`"bins"`** `int` Bins for projecting CVs at each dimension.
+
+Another example is given to project 2D path CV.
+```JSON
+{
+        "name": "rid-mcmc",
+        "init_models": ["model_000.pb","model_001.pb","model_002.pb","model_003.pb"],
+        
+        "MCMC_Config": {
+            "cv_dimension": 13,
+            "numb_steps": 100000,
+            "numb_walkers": 2000,
+            "cv_upper_bound": [5.33507,4.43402,4.59772,3.40287,4.67268,5.12813,3.28143,3.74967,3.58144,5.34143,7.13014,4.58024,7.41678],
+            "cv_lower_bound": [2.38722,0.99005,3.99886,2.08580,3.81799,3.93331,2.768175,3.145618,3.161886,4.661101,5.286805,4.076754,2.893047],
+            "proj_info": {
+                "proj_mode": "path",
+                "proj_cv_index": [0,1,2,3],
+                "path_list": [[2.38722,0.99005,3.99886,2.08580],[2.64067,1.20917,3.73831,1.86004],[2.90639,1.50430,3.47935,1.65185],
+                [3.18130,1.83919,3.22240,1.46871],[3.46322,2.19574,2.96795,1.32108],[3.75057,2.56491,2.71672,1.22191],
+                [4.04218,2.94198,2.46968,1.18344],[4.33720,3.32424,2.22824,1.21147],[4.63498,3.71010,1.99443,1.30170]],
+                "path_lm": 0.25
+            },
+            "cv_type": "dis",
+            "bins": 101
+        }
+    
+}
+```
+In the path CV projection mode, the additonal parameters are:
+* **`"proj_cv_index"`** `List[int]` Set of CV index used to defined the path CV.
+* **`"path_list"`** `List[float]` List of CVs to defined the path.
+* **`"path_lm"`** `float` Lambda in the path CV definition.
\ No newline at end of file
diff --git a/docs/source/prepare_files.md b/docs/source/prepare_files.md
new file mode 100644
index 00000000..f7ae876a
--- /dev/null
+++ b/docs/source/prepare_files.md
@@ -0,0 +1,29 @@
+# Preparing files for input to `RiD-kit`
+Currently Rid-kit manage files by their extensions, so you have to provide files with the correct extensions in the input directory. These includes the following:
+## For `Gromacs` simulation
+
+`.gro files` represent the initial configurations used by Gromacs run.
+
+`.top file` represent the topology information used by Gromacs run.
+
+`.ff directory (optional)` represents the forcefield information used by Gromacs, all the forcefield and restraint information not implemented in Gromacs natively should be included in this directory. Otherwise you can manually specify the forcefield information in the topology file.
+
+`.ndx file (optional)` represents the index file used by Gromacs, in situations where we want the control of different parts of the system.
+
+`.pb files (optional)` represents either the free energy model produced by Rid run, or the Deep Potential model used to run Gromacs. Specify different models in `rid.json` according to [Rid configuration](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/rid_configuration.md).
+
+`CV file (any extension or no extension) and its related .pdb files (optional)`: The CV file is used to customize CV definition, other files related to CV definition should has `pdb` extensions. The rules for defining customized CV can be found in [Rid configuration](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/rid_configuration.md).
+
+`.npy file (optional)` represents the CV-mean forces data used to train free energy model, this file is used when you want to use data from a previous RiD run.
+
+## For `Lammps` simulation
+
+`.lmp files` represent the initial configurations used by Lammps run.
+
+`input file (any extension or no extension)` represents the input file used by Lammps run.
+
+`.pb files (optional)` represents either the free energy model produced by Rid run, or the Deep Potential model used to run Gromacs. Specify different models in `rid.json` according to [Rid configuration](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/rid_configuration.md).
+
+`CV file (any extension or no extension) and its related .pdb files (optional)`: The CV file is used to customize CV definition, other files related to CV definition should has `pdb` extensions. The rules for defining customized CV can be found in [Rid configuration](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/rid_configuration.md).
+
+`.npy file (optional)` represents the CV-mean forces data used to train free energy model, this file is used when you want to use data from a previous RiD run.
\ No newline at end of file
diff --git a/docs/source/rid_configuration.md b/docs/source/rid_configuration.md
index fba57cd1..bddec5c0 100644
--- a/docs/source/rid_configuration.md
+++ b/docs/source/rid_configuration.md
@@ -205,12 +205,19 @@ This section configures the parameters in `Selection` step. In `Selection` step,
 This section configures the parameters in `Label` step. Currently, rid-kit supports two methods of labeling: `"restrained"` and `"constrained"`, stand for `restrained MD` and `constrained MD` respectively. Most settings are quite similar to those in `Exploration` Step. The simulation time is usually different between `Exploration` and `Label` steps, while the simulation time for `Exploration` has more freedom, the simulation time for `Label` step has to be chosen with care (longer enough to ensure convergence and shorter enough to avoid wasting). Our experience is that `100ps` is enough for `torsion mode` in `restrained MD method`, `1ns` is enough for `distance mode` in `constrained MD method`.  
 
 ### restrained method
-Set `"method": "restrained"` to use restrained MD as mean force calculator. The only different parameters with `Exploration` step is `kappas`.
+Set `"method": "restrained"` to use restrained MD as mean force calculator. The only different parameters with `Exploration` step is `kappas` and `std_threshold`.
 
 * **`kappas`** `(List[int])` A list of force constants ($\kappa$) of harmonic restraints. The length of the list is equal to the number of CVs.
 
+* **`std_threshold`** `(float)`(default 5.0, the unit is consistent with the mean force) A number represents the mean force standard deviation threshold, beyond which the mean force is neglected and will not be used in the dataset for training free energy model. You should test labeling MD for your own system to determine an appropriate number for this threshold.
+
 ### constrained method
-Set `"method": "constrained"` to use constrained MD as mean force calculator. Currently rid-kit only supports distance CV to use this method. The other parameters is the same with `Exploration` step.
+Set `"method": "constrained"` to use constrained MD as mean force calculator. Currently rid-kit only supports distance CV to use this method, also only `gmx` type is supported to perform constrained MD simulation. The other parameters is the same with `Exploration` step.
+
+Note that if you want to use constrained MD as the mean force calculator, apart from setting `method` to be `constrained` in the `label_config`, you should add `[ constraints ]` line corresponding to the `[ moleculartype ]` in your input `topology` file yourself, since gromacs specifies constraints information for each `[ moleculartype ]`.
+
+Also note that, since gromacs only supports constrained MD for `distance CV`, the constrained MD simulation in rid-kit only supports `distance CV` at this moment.
+
 
 ### Example
 
@@ -236,7 +243,8 @@ Set `"method": "constrained"` to use constrained MD as mean force calculator. Cu
         "ntmpi": 1,
         "nt": 8,
         "max_warning": 2,
-        "kappas": [ 500, 500 ]
+        "kappas": [ 500, 500 ],
+        "std_threshold": 2.0
     }
 ```
 ```JSON
@@ -260,7 +268,8 @@ Set `"method": "constrained"` to use constrained MD as mean force calculator. Cu
         "output_mode": "both",
         "ntmpi": 1,
         "nt": 8,
-        "max_warning": 2
+        "max_warning": 2,
+        "std_threshold": 10.0
 }
 ```
 
@@ -304,95 +313,6 @@ This section configures the parameters in `Train` step. `RiD-kit` is based on `T
 }
 ```
 
-## A full Example of `rid.json`
-
-You can find a full example of `rid.json` within `"rid-kit/rid/template"`. Or you can copy one from following:
+## Full Examples of `rid.json`
 
-```JSON
-{
-    "name": "test",
-    "numb_walkers": 2,
-    "numb_iters": 2,
-    "trust_lvl_1": 2,
-    "trust_lvl_2": 3,
-    "init_models": [],
-    
-    "CV": {
-        "mode": "torsion",
-        "selected_resid": [ 1, 2 ],
-        "angular_mask": [ 1, 1 ],
-        "weights": [ 1, 1 ],
-        "cv_file":[""]
-    },
-
-    "ExploreMDConfig": {
-        "nsteps": 50000,
-        "type": "gmx",
-        "temperature": 300,
-        "output_freq": 500,
-        "ref-t": "300 300",
-        "verlet-buffer-tolerance":"-1",
-        "rlist": 1,
-        "rvdw": 0.9,
-        "rvdw-switch": 0,
-        "rcoulomb": 0.9,
-        "rcoulomb-switch": 0,
-        "epsilon-r":1,
-        "epsilon-rf":80,
-        "dt": 0.002,
-        "fourierspacing": "0.12",
-        "output_mode": "single",
-        "ntmpi": 1,
-        "nt": 8,
-        "max_warning": 2
-    },
-
-    "SelectorConfig": {
-        "numb_cluster_lower": 12,
-        "numb_cluster_upper": 22,
-        "cluster_threshold": 1.5,
-        "max_selection": 30,
-        "numb_cluster_threshold": 6,
-        "slice_mode": "gmx"
-    },
-
-    "LabelMDConfig": {
-        "nsteps": 50000,
-        "temperature":300,
-        "method": "restrained",
-        "type": "gmx",
-        "output_freq": 100,
-        "ref-t": "300 300",
-        "rlist": 1,
-        "verlet-buffer-tolerance":"-1",
-        "rvdw": 0.9,
-        "rvdw-switch": 0,
-        "rcoulomb": 0.9,
-        "rcoulomb-switch": 0,
-        "epsilon-r":1,
-        "epsilon-rf":80,
-        "dt": 0.002,
-        "fourierspacing": "0.12",
-        "output_mode": "single",
-        "ntmpi": 1,
-        "nt": 8,
-        "max_warning": 2,
-        "kappas": [ 500, 500 ]
-    },
-
-    "Train": {
-        "numb_models": 4,
-        "neurons": [ 200, 200, 200, 200 ],
-        "resnet": true,
-        "batch_size": 128,
-        "epoches": 8000,
-        "init_lr": 0.0008,
-        "decay_steps": 120,
-        "decay_rate": 0.96,
-        "drop_out_rate": 0,
-        "numb_threads": 8,
-        "use_mix": false,
-        "restart": false
-    }
-}
-```
\ No newline at end of file
+You can find full examples of `rid.json` within `"rid-kit/rid/template"`.
\ No newline at end of file
diff --git a/docs/source/rid_download.md b/docs/source/rid_download.md
new file mode 100644
index 00000000..65a55a1d
--- /dev/null
+++ b/docs/source/rid_download.md
@@ -0,0 +1,46 @@
+# Downloading files in `RiD-kit`
+
+The download command is executed through
+```bash
+rid download WORKFLOW_ID -p STEP-KEY -f FILE_NAME -a ITERATION_START -e ITERATION_END -o OUTPUT_DIR
+```
+
+`ITERATION_START`(default 1) specifies from which iteration to start download, `ITERATION_END`(default 100) specifies to which iteration to end download. `OUTPUT_DIR`(default "./") is the output directories of the download.
+
+The `STEP-KEY` in rid includes the following steps: `prep-exploration`, `run-exploration`, `prep-select`, `run-select`, `prep-label`, `run-label`, `label-stats`, `collect-data`, `merge-data`, `train`, `model-devi`.
+
+The `FILE_NAME` in `prep-exploration` includes `task_path`. Which represents the input directory for Exploration MD simulation.
+
+The `FILE_NAME` in `run-exploration` includes `trajectory`, `conf_out`,`plm_out`,`md_log`,`bias_fig`, `model_devi_fig`, `dp_model_devi_fig (optional)`, `dp_model_devi (optional)`, `dp_selected_indices (optional)`, `dp_selected_confs (optional)`, `projected_fig (optional)`.
+
+The `trajectory` represents the trajectories of each walker produced by the Exploration MD, `conf_out` represents the last frame in the trajectories and will be used in the next iteration, `plm_out` represents the CV output from PLUMED, `md_log` is the log file from the molecular dynamics simulation, `bias_fig` is the figure showing the biasing potential used in the simulation, `model_devi_fig` is the figure showing the deviation of the free energy model used in the simulation, `dp_model_devi_fig (optional)` is the figure showing the deviation of the Deep Potential (DP) model used in the simulation, `dp_model_devi (optional)` contains numerical data on the deviation of the Deep Potential model used in the simulation, `dp_selected_indices (optional)` contains the indices of configurations with high dp model deviation selected for further analysis or training, `dp_selected_confs (optional)` contains the actual configurations high dp model deviation, `projected_fig (optional)` is the figure showing some projection of the simulation trajectory, currently only support for projection on the path CV. 
+
+The `FILE_NAME` in `prep-select` includes `cluster_fig`, `cluster_selection_index`, `cluster_selection_data`.
+
+The `cluster_fig`  is the figure showing the distribution of the clusters on the trajectories, `cluster_selection_index` contains the indices of the configurations from each cluster, `cluster_selection_data` contains the actual configurations selected from each cluster.
+
+The `FILE_NAME` in `run-select` includes `selected_confs`, `selected_cv_init`, `model_devi`, `selected_indices`,`selected_conf_tags`.
+
+The `selected_confs` contains the configurations that have been selected for labeling MD in the next step, `selected_cv_init` contains the initial values of the collective variables (CVs) for the selected configurations, `model_devi` contains the deviation of the free energy model on the selected configurations from the last `prep-select` step, `selected_indices` contains the indices of the selected configurations, `selected_conf_tags`  contains tags for the selected configurations. 
+
+The `FILE_NAME` in `prep-label` includes `task_path`. Which represents the input directory for Exploration MD simulation.
+
+The `FILE_NAME` in `run-label` includes `plm_out`, `cv_forces`, `mf_info`, `mf_fig`,`md_log`.
+
+The `plm_out` represents the CV output from PLUMED, `cv_forces` contains the collective variables (CVs) and the mean forces obtained from the simulation, `mf_info` contains information about the mean force (MF) calculated during the simulation, including the values and the standard deviation, `mf_fig` is the figure showing the running average of the mean force as a function of time, `md_log`  is a log file from the MD simulation.
+
+The `FILE_NAME` in `label-stats` includes `mf_std_fig` and `cv_forces`.
+
+The `mf_std_fig` is a histogram showing the distribution of the standard deviation of the mean force, `cv_forces` is a filtered version of the same file in last `run-label` step, removing mean forces with high standard deviation.
+
+The `FILE_NAME` in `collect-data` includes `data_new`. Which is the .npy file combining all the `cv_forces` files in the last step, represents the new training data.
+
+The `FILE_NAME` in `merge-data` includes `data_raw`. Which is the .npy file containing all the training data in the previous iterations.
+
+The `FILE_NAME` in `train` includes `model`, `train_log`, `train_fig`.
+
+The `model` is the trained free energy model, `train_log` is the log of the training process. `train_fig` represents the training error during the training process.
+
+The `FILE_NAME` in `model_devi` includes `model_devi`, `model_devi_fig`.
+
+The `model_devi` is the model deviation on the trajectories of each walker using the newly trained free energy model, `model_devi_fig` is the figure of the model deviation.
\ No newline at end of file
diff --git a/docs/source/troubleshooting/installation.md b/docs/source/troubleshooting/installation.md
index becad95b..e08cee62 100644
--- a/docs/source/troubleshooting/installation.md
+++ b/docs/source/troubleshooting/installation.md
@@ -10,6 +10,16 @@ add-apt-repository -y ppa:ubuntu-toolchain-r/test
 apt-get update
 apt-get install -y libstdc++6
 ```
+Or in slurm machine you might encounter this error in conda enviroment
+```
+/data1/ddwang/software/anaconda3/envs/rid-drop/lib/libgrpc.so.20: undefined reference to `std::__throw_bad_array_new_length()@GLIBCXX_3.4.29'
+/data1/ddwang/software/anaconda3/envs/rid-drop/lib/libtensorflow_framework.so.2: undefined reference to `std::__cxx11::basic_ostringstream<char, std::char_traits<char>, std::allocator<char> >::basic_ostringstream()@GLIBCXX_3.4.26'
+/data1/ddwang/software/anaconda3/envs/rid-drop/lib/libtensorflow_cc.so.2: undefined reference to `std::__cxx11::basic_stringstream<char, std::char_traits<char>, std::allocator<char> >::basic_stringstream()@GLIBCXX_3.4.26'
+```
+Try to install libgcc in conda
+```
+conda install libgcc
+```
 
 **You may encounter the error when compiling gromacs**
 ```
diff --git a/docs/source/tutorial.ipynb b/docs/source/tutorial.ipynb
index 8e164da9..fa9cf8c3 100644
--- a/docs/source/tutorial.ipynb
+++ b/docs/source/tutorial.ipynb
@@ -1,6 +1,7 @@
 {
  "cells": [
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -10,6 +11,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -22,6 +24,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -29,6 +32,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -36,6 +40,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -43,6 +48,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -72,6 +78,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -89,6 +96,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -96,6 +104,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -107,6 +116,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -126,6 +136,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -152,6 +163,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -181,6 +193,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -229,6 +242,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -249,6 +263,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -264,10 +279,44 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "! rid submit -i /tests/data/000 -c /rid/template/rid_gmx_dih.json -m /rid/template/machine_bohrium_local.json"
+    "! rid submit -i /tests/data/000 -c /rid/template/rid_gmx_dih.json -m /rid/template/machine_bohrium_k8s.json"
    ]
   },
   {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can specify the workflow name by providing WORKFLOW_ID after \"-d\", for example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "! rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_bohrium_k8s.json -d ala-dipeptide-1"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Note that the defined workflow-id should only contain lower case alphanumeric character, and specifal character \"-\".**"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that the defined workflow-id should only contain lower case alphanumeric character, and specifal character \"-\"."
+   ]
+  },
+  {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -277,6 +326,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -293,6 +343,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -300,6 +351,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -319,6 +371,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -331,10 +384,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "! rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json Workflow-ID"
+    "# suppose the original workflow id is OLD_ID\n",
+    "! rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json OLD_ID -d NEW_ID"
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -347,14 +402,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "! rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json Workflow-ID -t ITERATION-ID -p STEP-KEY"
+    "! rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json OLD_ID -t ITERATION-ID -p STEP-KEY -d NEW_ID"
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The `Workflow-ID` is something like `reinforced-dynamics-jq4jn` appeared in the argo UI. If the workflow is archived, its name will appear as somethinig like `a8463748-e15a-4f2c-882b-bfd981a76dac`. If you want to resubmit a archived workflow, you have to provide the archived ID rather than the initial ID. The `ITERATION-ID` is just `n`th iteration the workflow has been executed. The `STEP-KEY` in rid includes the following steps: `prep-exploration`, `run-exploration`, `prep-select`, `run-select`, `prep-label`, `run-label`, `post-label`, `collect`, `merge`, `train`."
+    "The `Workflow-ID` is something like `reinforced-dynamics-jq4jn` appeared in the argo UI. If the workflow is archived, its name will appear as somethinig like `a8463748-e15a-4f2c-882b-bfd981a76dac`. If you want to resubmit a archived workflow, you have to provide the archived ID rather than the initial ID. The `ITERATION-ID` is just `n`th iteration the workflow has been executed. The `STEP-KEY` in rid includes the following steps: `prep-exploration`, `run-exploration`, `prep-select`, `run-select`, `prep-label`, `run-label`, `label-stats`, `collect-data`, `merge-data`, `train`, `model-devi`."
    ]
   }
  ],
diff --git a/docs/source/tutorial.md b/docs/source/tutorial.md
new file mode 100644
index 00000000..919594ed
--- /dev/null
+++ b/docs/source/tutorial.md
@@ -0,0 +1,195 @@
+# Tutorial
+
+In this tutorial, we will learn how to deploy `RiD-kit` with your own `Kubenetes` enviroment and run a simple case of alanine dipeptide.
+
+
+## Installation of `dflow` and `rid-kit`
+With the power of `dflow`, users can easily minitor the whole workflow of RiD tasks and dispatch their tasks to various computational resources. Before you use it, you should have `dflow` installed on your host computer (your PC or a remote server). 
+
+It it necessary to emphasize that, the computational nodes and monitor nodes are seperated. With `dflow`, you can deploy `dflow` and `rid` on your PC and achieve expensive computation on other resources (like `Slurm` and Cloud Platform) without any further effort.
+
+Instructions of `dflow` installation are provided in detail on its [Github page](https://github.com/deepmodeling/dflow#Installdflow). Prerequisites of `dflow` usage are `Docker` and `Kubenetes`, where their main pages ([Docker](https://docs.docker.com/engine/install/) & [Kubenetes](https://kubernetes.io/docs/tasks/tools/)) include how you can install them. Besides, `dflow` repo also provides with easy-install shell scripts on [dflow/scripts](https://github.com/deepmodeling/dflow/tree/master/scripts) to install `Docker` & `Kubenetes` & `dflow` and make port-forwarding.
+
+Instead of using the easy-install scripts, it is recommended to configure your own minikube enviroment
+
+Installation of `minikube` is very easy, just follow the page [Minikube](https://minikube.sigs.k8s.io/docs/start/) based on your machine archetecture.
+
+Then use minikube to start a k8s environment. This can be done by just type `minikube start`, but since running rid-kit requires relatively huge amount of resources, it is recommended to allocate more resources when executing `minikube start`, just as follows.
+
+> **Note:**
+>
+> Don't try to run `minikube` with root privileges, otherwise an error may occur:
+>
+> `Exiting due to DRV_AS_ROOT: The "docker" driver should not be used with root privileges.`
+>
+
+```bash
+#choose the location of the minikube enviroment
+export MINIKUBE_HOME=~/.minikube  
+#allocate enough memory based on your machine in case of high parallelism
+minikube start --cpus 8 --memory 8192mb --kubernetes-version=1.23.9 --image-mirror-country='cn'
+# mount the storage path if you are using machine with shared memory storage system, change the minio host path accordingly
+minikube start --cpus 8 --memory 8192mb --kubernetes-version=1.23.9 --mount --mount-string="/path_on_your_machine:/data2" --image-mirror-country='cn'
+```
+
+A further step to  configure `argo` service is to run:
+
+```bash
+kubectl create ns argo
+kubectl apply -n argo -f https://raw.githubusercontent.com/deepmodeling/dflow/master/manifests/quick-start-postgres.yaml
+```
+
+If you want to mount your local path into the minio path, change the minio host path in the yaml file
+
+```yaml
+      - hostPath:
+          path: /data2/minio
+          type: DirectoryOrCreate
+```
+
+Now you should have installed `Docker` and `minikube` properly. Run commands to check their status. For `minikube`, you should wait util all servers keep `running`. This may take a couple of minutes.
+
+```bash
+# check minikube status
+minikube status
+# check k8s status in argo namespace
+kubectl get pods -n argo
+```
+
+### Installation of RiD-kit
+
+Now we install `rid-kit` on the host machine. To meet the minimum requirments, the prerequisites of third-party python package should be installed:
+
+* tensorflow-cpu or gpu
+* mdtraj
+* numpy
+* scikit-learn
+* pydflow
+* google
+* mdtraj
+* dpdata
+* parmed
+* dpdispatcher
+* lbg
+* matplotlib
+
+which are also listed in `rid-kit/requirements.txt`. Then change directory to `rid-kit` repo and run:
+
+```bash
+# the rid-kit repo path
+git checkout dflow
+pip install setuptools_scm
+pip install .
+```
+
+## Configuration of Computational Environment
+
+In RiD workflow, `dflow` helps send computation tasks to resources with peoper environment configured. 
+
+There are four main modules and several workflow steps in RiD procedures and each module or step needs different environments:
+
+* `Exploration/Sampling`:  `Gromacs`, `PLUMED2` modified by `DeepFE.cpp`, `Tensorflow` C++ interface. (prefer GPU)
+* `Selection`:  `Tensroflow` Python interface.
+* `Labeling`:  `Gromacs`, `PLUMED2`. (prefer GPU)
+* `Training`:  `Tensroflow` Python interface. (prefer GPU).
+* `Workflow steps`:  `Python`.
+
+`dflow` supports different resources including `Slurm` clusters, `K8S` local machines and `Cloud Server`.
+
+* For `Slurm`, configure computational environments on your `Slurm` following the instructions in [Environment settings](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/install.md). With `dflow`, `rid-kit` send tasks to `Slurm` nodes from the host machines remotely without manually logging in the cluster.
+* For local resources, just use the docker images we have built. No further manual configuration needed. We also provide `Dockerfile` of our images to enable flexible modification.
+* For `Cloud Server`, like `Bohrium`, use public images and no further manual configuration needed.
+
+We highly recommend using `Bohrium` to do computation since no enviroment installation is required, you only need to use the public image.Specifically, the public images for rid-kit involves the following: `registry.dp.tech/public/pkufjhdocker/rid-gmx-exploration:stable` for doing `Exploration`, `registry.dp.tech/public/pkufjhdocker/rid-gmx-plumed:stable` for doing `Labeling`, `registry.dp.tech/public/pkufjhdocker/rid-gmx-tf:stable` for doing `Selection`, `registry.dp.tech/public/pkufjhdocker/rid-tf-gpu:stable` for `Training` and `registry.dp.tech/public/pkufjhdocker/rid-tf-cpu:stable` for doing other steps.
+
+## Prepare machine configuration JSON.
+
+`rid-kit` uses `JSON` file to manage resources. In `machine.json`, define your own `resources` and dispatch `tasks` to them.
+* In key `resources`, you define your own resources types. Resource names and their numbers are custom.
+* In key `tasks`, you distribute resources you have defined to tasks of RiD. Do not change task names in `tasks` as they are fixed in codes.
+
+Generally, we would like to run low-cost tasks on cpu nodes or locally and submit high-cost tasks to `Slurm` or `Clouds`. 
+
+If submitting to mixed `Slurm` and `local` environment, a `machine.json` example is in [rid/template/machine_slurm_k8s.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/machine_slurm_k8s.json):
+
+If submitting to `Bohrium` cloud platform, a `machine.json` example is in [rid/template/machine_bohrium.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/machine_bohrium.json).
+
+If submitting to mixed `Bohrium` and `local` platform, a `machine.json` example is in [rid/template/machine_bohrium_k8s.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/machine_bohrium_k8s.json).
+
+If submitting to `Slurm` environment without `k8s` enviroment, a `machine.json` example is in [rid/template/machine_slurm_local.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/machine_slurm_local.json):
+
+If you submit jobs to `Slurm` enviroment, you will have to compile the computation enviroment on `Slurm` server for yourself. While in `Bohrium` cloud platform the computation images are all set, no extra compilation is needed. So we recommend to submit jobs to `Cloud` enviroment like `Bohrium`.
+
+For more details about setting the machine resources, check [Rid machine resourse](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/rid_machine.md)
+
+## Prepare rid configuration JSON.
+If you want to use `gromacs` as sampler, [rid/template/rid_gmx_dih.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/rid_gmx_dih.json) is an example.
+
+If you want to use `gromacs` with `DeepMD` potential as sampler, [rid/template/rid_gmx_dp.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/rid_gmx_dp.json) is an example.
+
+If you want to use `lammps` with `DeepMD` potential as sampler, [rid/template/rid_lmp_dp.json](https://github.com/PKUfjh/rid-kit/blob/dflow/rid/template/rid_lmp_dp.json) is an example.
+
+For more details about setting the rid configuration, check [Rid configuration](https://github.com/PKUfjh/rid-kit/blob/dflow/docs/source/rid_configuration.md)
+
+## Get Started!
+
+Assume you have learn the basic knowledge of reinforced dynamics which we won't describe again here. 
+
+Users can monitor workflows from browser UI. To enable that, you should forward ports of `argo` and `minio`. These could be achieved by `rid port-forward`.
+
+```bash
+rid port-forward
+```
+
+In this case, we try to explore the phase space of alanine dipeptide. Prepare your initial conformation files in `.gro` format, topology file in `.top` format and configuration file `rid.json`. For convenience, we have prepared on at `rid/template/rid.json`. 
+Remember also provide your own forcefield files. Collect all these files into a directory and feed its path to `rid-kit` by flag `-i`.
+
+A minimum case was prepared in `rid-kit/tests/data/000`. After configurating your `machine.json`, then run `rid submit`:
+
+```bash
+rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_bohrium_k8s.json
+```
+
+You can specify the workflow name by providing WORKFLOW_ID after "-d", for example:
+
+```bash
+rid submit -i ./tests/data/000 -c ./rid/template/rid_gmx_dih.json -m ./rid/template/machine_bohrium_k8s.json -d ala-dipeptide-1
+```
+
+**Note that the defined workflow-id should only contain lower case alphanumeric character, and specifal character "-".**
+
+INFO indicates that this task has been submitted succussfully. Record this workflow ID as we may use it later.
+
+Visit the `url` given by the last two lines, all workflows and corresponding files are listed on UI.
+
+Command lines are also supported. Run `rid ls` to list your workflows and their status.
+
+```bash
+! rid ls
+```
+
+`rid-kit` is based on `dflow`, `argo` and `minikube`. So further complex and flexible managements of workflows can be achieved by their command lines. like `kubectl get pods -n argo` and `argo show`.
+
+For failed tasks, you may want to remove them or resubmit them from the failure steps.
+
+For `remove`:
+
+```bash
+# rid rm task-ID
+rid rm reinforced-dynamics-bsc7j 
+```
+
+For `resubmit` to modify and continue workflow:
+
+```bash
+# suppose the original workflow id is OLD_ID
+rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json OLD_ID -d NEW_ID
+```
+
+If you want to resubmit from a particular `iteration` and `step`:
+
+```bash
+rid resubmit -i your_dir -c path_to_rid.json -m path_to_machine.json OLD_ID -t ITERATION-ID -p STEP-KEY -d NEW_ID
+```
+
+The `Workflow-ID` is something like `reinforced-dynamics-jq4jn` appeared in the argo UI. If the workflow is archived, its name will appear as somethinig like `a8463748-e15a-4f2c-882b-bfd981a76dac`. If you want to resubmit a archived workflow, you have to provide the archived ID rather than the initial ID. The `ITERATION-ID` is just `n`th iteration the workflow has been executed. The `STEP-KEY` in rid includes the following steps: `prep-exploration`, `run-exploration`, `prep-select`, `run-select`, `prep-label`, `run-label`, `label-stats`, `collect-data`, `merge-data`, `train`, `model-devi`.
diff --git a/install/lmp_exploration_dp.dockerfile b/install/lmp_exploration_dp.dockerfile
index 70393ca2..d51892a2 100644
--- a/install/lmp_exploration_dp.dockerfile
+++ b/install/lmp_exploration_dp.dockerfile
@@ -30,7 +30,7 @@ RUN conda init bash && source /root/.bashrc \
     && conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/  \
     && conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/pytorch/  \
     && conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge  \
-    && conda install python=3.9 libtensorflow_cc=2.6.2=*cuda110* dpdata nccl numpy -y
+    && conda install python=3.9 libtensorflow_cc=2.6.2=*cuda110* tensorflow=2.6.2=*cuda110* dpdata nccl numpy -y
 
 # Solve library inconsistency
 RUN source activate base && rm ${CONDA_PREFIX}/lib/libtinfo.so* && ln -s /usr/lib/x86_64-linux-gnu/libtinfo.so.6 ${CONDA_PREFIX}/lib/libtinfo.so.6 \
diff --git a/requirements.txt b/requirements.txt
index 1781e309..c1d0b3f2 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,7 +1,8 @@
 protobuf==3.20.1
 scikit-learn
 pydflow
-tensorflow==2.8.0
+tensorflow==2.6.2
+pandas==1.4.4
 google
 mdtraj
 dpdata
diff --git a/rid/__init__.py b/rid/__init__.py
index acc3d8ba..405deae5 100644
--- a/rid/__init__.py
+++ b/rid/__init__.py
@@ -11,4 +11,5 @@
 SRC_ROOT = __path__[0]
 config["extender_image_pull_policy"] = "IfNotPresent"
 config["util_image_pull_policy"] = "IfNotPresent"
-config["dispatcher_image_pull_policy"] = "IfNotPresent"
\ No newline at end of file
+config["dispatcher_image_pull_policy"] = "IfNotPresent"
+config["save_keys_in_global_outputs"] = False
\ No newline at end of file
diff --git a/rid/common/mol.py b/rid/common/mol.py
index bc9f070c..c8152ad8 100644
--- a/rid/common/mol.py
+++ b/rid/common/mol.py
@@ -3,7 +3,6 @@
 from typing import List, Dict, Sequence, Union
 import logging
 import mdtraj as md
-import dpdata
 from mdtraj.geometry.dihedral import _atom_sequence, PHI_ATOMS, PSI_ATOMS
 import numpy as np
 from rid.constants import sel_gro_name_gmx, sel_gro_name
@@ -122,32 +121,5 @@ def slice_xtc(
             selected_idx = selected_idx,
             output_format=output
         )
-    else:
-        raise RuntimeError("Unknown Style for Slicing Trajectory.")
-
-def slice_dump_dpdata(
-    dump: str,
-    walker_idx,
-    selected_idx,
-    output_format:str
-):
-    traj = dpdata.System(dump, fmt="lammps/dump")
-    for sel in selected_idx:
-        traj.to('lammps/lmp', output_format.format(walker=walker_idx,idx=sel), frame_idx=sel)
-
-def slice_dump(
-    dump: str,
-    walker_idx: int,
-    selected_idx,
-    output:str,
-    style: str = "dpdata"
-):
-    if style == "dpdata":
-        slice_dump_dpdata(
-            dump = dump,
-            walker_idx= walker_idx,
-            selected_idx = selected_idx,
-            output_format=output
-        )
     else:
         raise RuntimeError("Unknown Style for Slicing Trajectory.")
\ No newline at end of file
diff --git a/rid/common/mol_dpdata.py b/rid/common/mol_dpdata.py
new file mode 100644
index 00000000..393315de
--- /dev/null
+++ b/rid/common/mol_dpdata.py
@@ -0,0 +1,41 @@
+import dpdata
+from typing import List, Dict, Sequence, Union
+
+def slice_dump_dpdata(
+    dump: str,
+    walker_idx,
+    selected_idx,
+    output_format:str,
+    type_map: List = [],
+):
+    if type_map != []:
+        traj = dpdata.System(dump, fmt="lammps/dump", type_map=type_map)
+    else:
+        traj = dpdata.System(dump, fmt="lammps/dump")
+    if output_format.endswith("lmp"):
+        for sel in selected_idx:
+            traj.to('lammps/lmp', output_format.format(walker=walker_idx,idx=sel), frame_idx=sel)
+    elif output_format.endswith("gro"):
+        for sel in selected_idx:
+            traj.to('gromacs/gro', output_format.format(walker=walker_idx,idx=sel), frame_idx=sel)
+    else:
+        raise ValueError("Invalid output format, only support gmx and lmp")
+
+def slice_dump(
+    dump: str,
+    walker_idx: int,
+    selected_idx,
+    output:str,
+    type_map: List = [],
+    style: str = "dpdata"
+):
+    if style == "dpdata":
+        slice_dump_dpdata(
+            dump = dump,
+            walker_idx= walker_idx,
+            selected_idx = selected_idx,
+            type_map = type_map,
+            output_format=output
+        )
+    else:
+        raise RuntimeError("Unknown Style for Slicing Trajectory.")
\ No newline at end of file
diff --git a/rid/constants.py b/rid/constants.py
index cf658bec..bffd99ad 100644
--- a/rid/constants.py
+++ b/rid/constants.py
@@ -1,13 +1,14 @@
 # rid task name format
 argo_namespace = "argo"
 walker_tag_fmt = "{idx:03d}"
-init_conf_name = "conf_{idx:03d}.gro"
+init_conf_gmx_name = "conf_{idx:03d}.gro"
+init_conf_lmp_name = "conf_{idx:03d}.lmp"
 init_input_name = "input_{idx:03d}.lammps"
 explore_task_pattern = "{:03d}"
 explore_task_file = "explore_{walker:03d}.pkl"
 cluster_selection_data_name = "cls_sel.out.npy"
 cluster_selection_index_name = "cls_sel.ndx.npy"
-sel_ndx_name = "sel.ndx"
+sel_ndx_name = "sel.ndx.npy"
 cv_init_label = "cv_init_{walker:03d}_{idx:d}.out"
 model_devi_name = "model_devi.txt"
 label_task_pattern = "{:03d}"
@@ -58,20 +59,33 @@
 
 # Dp config file
 dp_config_name = "dp_config"
+dp_sel_ndx = "dp_sel.ndx"
 
 # figure file name
 mf_fig = "mf_average.png"
 bias_fig = "bias.png"
 model_devi_fig = "model_devi.png"
+new_model_devi_fig = "new_model_devi.png"
 cluster_fig = "cluster.png"
 train_fig = "train_{tag}.png"
+dp_model_devi_fig = "dp_model_devi.png"
+mf_std_fig = "mf_std.png"
+
+# MCMC file name
+mcmc_1cv_dir_name = "mcmc_1cv"
+mcmc_2cv_name = "mcmc_2cv.dat"
+mcmc_1cv_name = "mcmc_1cv_{tag}.dat"
+mcmc_2cv_fig = "mcmc_2cv.png"
+mcmc_1cv_fig = "mcmc_1cv_{tag}.png"
+mcmc_2cv_fig_separate = "mcmc_2cv_{tag}.png"
 
 # Units
-kb = 8.617333E-5
+kb = 8.617333E-5 # ev
 kbT = (8.617333E-5) * 300
 beta = 1.0 / kbT
-f_cvt = 96.485
+f_cvt = 96.485 # ev to kj/mol
 inverse_f_cvt = 1 / f_cvt
+kcal2kj = 4.184
 
 # precision
 model_devi_precision = "%.6e"
\ No newline at end of file
diff --git a/rid/entrypoint/download.py b/rid/entrypoint/download.py
new file mode 100644
index 00000000..c38e6e98
--- /dev/null
+++ b/rid/entrypoint/download.py
@@ -0,0 +1,53 @@
+from dflow import Workflow, download_artifact
+import numpy as np
+import os
+from datetime import datetime
+
+from dflow import config, s3_config
+from dflow.plugins import bohrium
+from dflow.plugins.bohrium import TiefblueClient
+
+
+def download_rid(workflow_id, pod, file, iteration_start, iteration_end, outputs):
+    wf = Workflow(id=workflow_id)
+    all_steps = wf.query_step()
+
+    iter_a = int(iteration_start)
+    iter_e = int(iteration_end)+1
+    print("start from iteration",iter_a)
+    print("end to iteration",iter_e-1)
+    for step in all_steps:
+        if step["type"] == "Pod":
+            for index in range(iter_a, iter_e):
+                # print(step["key"])
+                if step['key'] == "iter-%03d-%s-0"%(index,pod):
+                    print("downloading "+step["key"])
+                    download_artifact(step.outputs.artifacts[file], path="%s/iter%s"%(outputs,index))
+                elif step['key'] == "iter-%03d-%s"%(index,pod):
+                    print("downloading "+step["key"])
+                    download_artifact(step.outputs.artifacts[file], path="%s/iter%s"%(outputs,index))
+                elif step["key"] == "000-%s-0"%pod:
+                    print("downloading "+step["key"])
+                    download_artifact(step.outputs.artifacts[file], path="%s"%outputs)
+                    break
+                elif step["key"] == "000-%s"%pod:
+                    print("downloading "+step["key"])
+                    download_artifact(step.outputs.artifacts[file], path="%s"%outputs)
+                    break
+
+# with open("time_profile.txt","w") as f:
+#     for step in all_steps:
+#         if step['key'] is not None:
+#             f.write("step key:"+str(step["key"]))
+#         if step['startedAt'] is not None:
+#             start_time_str = str(step['startedAt'])
+#             f.write(" start time:"+str(start_time_str))
+#             start_time = datetime.strptime(start_time_str, "%Y-%m-%dT%H:%M:%SZ")
+#         if step['finishedAt'] is not None:
+#             end_time_str = str(step['finishedAt'])
+#             f.write(" finish time:"+str(end_time_str))
+#             end_time = datetime.strptime(end_time_str, "%Y-%m-%dT%H:%M:%SZ")
+#             # Calculate the duration of the task
+#             duration = end_time - start_time
+#             f.write(" duration:"+str(duration))
+#         f.write("\n")
\ No newline at end of file
diff --git a/rid/entrypoint/explore.py b/rid/entrypoint/explore.py
new file mode 100644
index 00000000..2863e574
--- /dev/null
+++ b/rid/entrypoint/explore.py
@@ -0,0 +1,196 @@
+from distutils.command.config import dump_file
+import json
+from pathlib import Path
+from typing import List, Union, Optional
+from rid.utils import load_json
+from copy import deepcopy
+import os
+
+from dflow import (
+    Workflow,
+    Step,
+    upload_artifact
+)
+
+from dflow.python import(
+    PythonOPTemplate,
+    upload_packages,
+    OP
+)
+from rid import SRC_ROOT
+upload_packages.append(SRC_ROOT)
+
+from rid.utils import normalize_resources,init_executor
+from rid.op.prep_rid import PrepRiD
+from rid.superop.exploration import Exploration
+from rid.op.prep_exploration import PrepExplore
+from rid.op.run_exploration import RunExplore
+
+def explore_rid(
+        confs: Union[str, List[str]],
+        topology: Optional[str],
+        rid_config: str,
+        machine_config: str,
+        models: Optional[Union[str, List[str]]] = None,
+        forcefield: Optional[str] = None,
+        index_file: Optional[str] = None,
+        dp_files: Optional[List[str]] = [],
+        otherfiles: Optional[List[str]] = None
+    ):
+    with open(machine_config, "r") as mcg:
+        machine_config_dict = json.load(mcg)
+    resources = machine_config_dict["resources"]
+    tasks = machine_config_dict["tasks"]
+    normalized_resources = {}
+    for resource_type in resources.keys():
+        normalized_resources[resource_type] = normalize_resources(resources[resource_type])
+        
+    if isinstance(confs, str):
+        confs_artifact = upload_artifact(Path(confs), archive=None)
+    elif isinstance(confs, List):
+        confs_artifact = upload_artifact([Path(p) for p in confs], archive=None)
+    else:
+        raise RuntimeError("Invalid type of `confs`.")
+
+    if index_file is None:
+        index_file_artifact = None
+    else:
+        index_file_artifact = upload_artifact(Path(index_file), archive=None)
+    
+    jdata = deepcopy(load_json(rid_config))
+    cv_config = jdata["CV"]
+    explore_config = jdata["ExploreMDConfig"]
+    
+    inputfiles = []
+    if "inputfile" in explore_config:
+        inputfiles.append(explore_config["inputfile"])
+
+    fe_models = []
+    assert isinstance(jdata["init_models"],list), "model input should be list."
+    for model in jdata["init_models"]:
+        fe_models.append(model)
+    
+    cvfiles = []
+    if "cv_file" in jdata["CV"]:
+        assert isinstance(jdata["CV"]["cv_file"],list), "CV file input should be list."
+        for file in jdata["CV"]["cv_file"]:
+            cvfiles.append(file)
+            
+    dp_models = []
+    if "dp_model" in explore_config:
+        assert isinstance(explore_config["dp_model"],list), "model input should be list."
+        for model in explore_config["dp_model"]:
+            dp_models.append(model)
+                
+    inputfile_list = []
+    cvfile_list = []
+    model_list = []
+    dpfile_list = []
+    if otherfiles is not None:
+        for file in otherfiles:
+            if os.path.basename(file) in inputfiles:
+                inputfile_list.append(file)
+            elif os.path.basename(file) in cvfiles:
+                cvfile_list.append(file)
+                
+    if dp_files is not None:
+        for dp_file in dp_files:
+            dpfile_list.append(dp_file)
+            
+    if models is not None:
+        for model in models:
+            if os.path.basename(model) in fe_models:
+                model_list.append(model)
+            elif os.path.basename(model) in dp_models:
+                dpfile_list.append(model)
+            
+    if len(inputfile_list) == 0:
+        inputfile_artifact = None
+    else:
+        inputfile_artifact = upload_artifact([Path(p) for p in inputfile_list], archive=None)
+        
+    if len(model_list) == 0:
+        models_artifact = None
+    else:
+        models_artifact = upload_artifact([Path(p) for p in model_list], archive=None)
+        
+    if len(cvfile_list) == 0:
+        cv_file_artifact = None
+    else:
+        cv_file_artifact = upload_artifact([Path(p) for p in cvfile_list], archive=None)
+        
+    if len(dpfile_list) == 0:
+        dp_files_artifact = None
+    elif isinstance(dp_files, List):
+        dp_files_artifact = upload_artifact([Path(p) for p in dpfile_list], archive=None)
+    else:
+        raise RuntimeError("Invalid type of `dp_files`.")
+    
+    if forcefield is None:
+        forcefield_artifact = None
+    else:
+        forcefield_artifact = upload_artifact(Path(forcefield), archive=None)
+        
+    if topology is None:
+        top_artifact = None
+    else:
+        top_artifact = upload_artifact(Path(topology), archive=None)
+        
+    rid_config = upload_artifact(Path(rid_config), archive=None)
+        
+    _step_config = normalized_resources[tasks["workflow_steps_config"]]
+    _step_config = deepcopy(_step_config)
+    step_template_config = _step_config.pop('template_config')
+    step_executor = init_executor(_step_config.pop('executor'))
+    
+    prep_rid = Step(
+        name = 'prepare-rid',
+        template=PythonOPTemplate(
+            PrepRiD,
+            python_packages = None,
+            **step_template_config,
+        ),
+        parameters={},
+        artifacts={
+            "confs":confs_artifact,
+            "rid_config" : rid_config
+        },
+        key = 'prepare-rid',
+        executor = step_executor,
+        **_step_config,
+    )
+
+    explore_op = Exploration(
+        "explore",
+        PrepExplore,
+        RunExplore,
+        prep_config = normalized_resources[tasks["prep_exploration_config"]],
+        run_config = normalized_resources[tasks["run_exploration_config"]],
+        retry_times=1)
+
+
+    exploration_steps = Step("rid-explore",
+            explore_op,
+            artifacts={
+                "models": models_artifact,
+                "forcefield": forcefield_artifact,
+                "topology": top_artifact,
+                "inputfile": inputfile_artifact,
+                "confs": prep_rid.outputs.artifacts["confs"],
+                "index_file": index_file_artifact,
+                "dp_files": dp_files_artifact,
+                "cv_file": cv_file_artifact
+            },
+            parameters={
+                "trust_lvl_1" : prep_rid.outputs.parameters["trust_lvl_1"],
+                "trust_lvl_2": prep_rid.outputs.parameters["trust_lvl_2"],
+                "exploration_config" :explore_config,
+                "cv_config": cv_config,
+                "task_names": prep_rid.outputs.parameters["walker_tags"],
+                "block_tag" : "iter-001"
+            },
+        )
+    wf = Workflow("rid-exploration", pod_gc_strategy="OnPodSuccess", parallelism=50)
+    wf.add(prep_rid)
+    wf.add(exploration_steps)
+    wf.submit()
diff --git a/rid/entrypoint/main.py b/rid/entrypoint/main.py
index fccb7cc2..f37bda0a 100644
--- a/rid/entrypoint/main.py
+++ b/rid/entrypoint/main.py
@@ -27,8 +27,13 @@
 
 from .submit import submit_rid
 from .resubmit import resubmit_rid
+from .explore import explore_rid
 from .label import label_rid
 from .relabel import relabel_rid
+from .redim import redim_rid
+from .reredim import reredim_rid
+from .train import train_rid
+from .download import download_rid
 from .info import information
 from .server import forward_ports
 from .cli import rid_ls, rid_rm
@@ -107,6 +112,9 @@ def main_parser() -> argparse.ArgumentParser:
     parser_run.add_argument(
         "--machine", "-m", help="Machine configuration.", dest="machine"
     )
+    parser_run.add_argument(
+        "--id", "-d", help="Workflow user defined ID", default = None, dest="id"
+    )
 
     # resubmit
     parser_rerun = subparsers.add_parser(
@@ -126,6 +134,9 @@ def main_parser() -> argparse.ArgumentParser:
     parser_rerun.add_argument(
         "--machine", "-m", help="Machine configuration.", dest="machine"
     )
+    parser_rerun.add_argument(
+        "--id", "-d", help="Workflow user defined ID", default = None, dest="id"
+    )
     parser_rerun.add_argument(
         "--iteration", "-t", help="restart from t-th iteration.", default = None, dest="iteration"
     )
@@ -135,8 +146,8 @@ def main_parser() -> argparse.ArgumentParser:
     
     # Explore
     parser_exp = subparsers.add_parser(
-        "mdrun",
-        help="Submit RiD workflow",
+        "explore",
+        help="Exploration MD",
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
     )
     parser_exp.add_argument(
@@ -191,10 +202,13 @@ def main_parser() -> argparse.ArgumentParser:
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
     )
     parser_train.add_argument(
-        "--data", "-d", help="Training data."
+        "--mol", "-i", help="Training data.", dest="mol",
+    )
+    parser_train.add_argument(
+       "--config", "-c", help="RiD configuration.", dest="config"
     )
     parser_train.add_argument(
-        "--config", "-c", help="RiD configuration."
+         "--machine", "-m", help="Machine configuration.", dest="machine"
     )
 
     # NN dimension reduction.
@@ -204,13 +218,65 @@ def main_parser() -> argparse.ArgumentParser:
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
     )
     parser_redim.add_argument(
-        "--networks", "-n", help="Training data."
+        "--mol", "-i", help="Neural networks path", dest="mol",
+    )
+    parser_redim.add_argument(
+        "--config", "-c", help="RiD configuration.", dest="config"
     )
     parser_redim.add_argument(
-        "--dim1", help="Dimension 1."
+        "--machine", "-m", help="Machine configuration.", dest="machine"
     )
     parser_redim.add_argument(
-        "--dim2", help="Dimension 2."
+        "--id", "-d", help="Workflow user defined ID", default = None, dest="id"
+    )
+    
+    # resubmit NN dimension reduction.
+    parser_reredim = subparsers.add_parser(
+        "reredim",
+        help="NN dimension reduction by Monte Carlo integral.",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    )
+    parser_reredim.add_argument(
+        "--mol", "-i", help="Neural networks path", dest="mol",
+    )
+    parser_reredim.add_argument(
+        "--config", "-c", help="RiD configuration.", dest="config"
+    )
+    parser_reredim.add_argument(
+        "--machine", "-m", help="Machine configuration.", dest="machine"
+    )
+    parser_reredim.add_argument(
+        "--id", "-d", help="Workflow user defined ID", default = None, dest="id"
+    )
+    parser_reredim.add_argument(
+        "WORKFLOW_ID", help="Workflow ID."
+    )
+    parser_reredim.add_argument(
+        "--pod", "-p", help="restart from the pod.", default = None, dest="pod"
+    )
+    # download
+    parser_download = subparsers.add_parser(
+        "download",
+        help="Download files in RiD workflow",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    )
+    parser_download.add_argument(
+        "WORKFLOW_ID", help="Workflow ID."
+    )
+    parser_download.add_argument(
+        "--pod", "-p", help="The pod in RiD workflow", dest="pod"
+    )
+    parser_download.add_argument(
+        "--file", "-f", help="The file to download inside the step", dest = "file"
+    )
+    parser_download.add_argument(
+        "--iteration_start", "-a", help="download from a-th iteration.", default = 1, dest="iteration_start"
+    )
+    parser_download.add_argument(
+        "--iteration_end", "-e", help="download to e-th iteration.", default = 100, dest="iteration_end"
+    )
+    parser_download.add_argument(
+        "--outputs", "-o", help="output directory of downloading", default = "./", dest="outputs"
     )
 
     # --version
@@ -249,6 +315,7 @@ def parse_submit(args):
     top_file = None
     forcefield = None
     index_file = None
+    plm_files = []
     dp_files = []
     models = []
     data_file = None
@@ -268,27 +335,33 @@ def parse_submit(args):
             index_file = file
         elif os.path.basename(file).endswith("json") or os.path.basename(file).endswith("raw"):
             dp_files.append(file)
+        elif os.path.basename(file).endswith("out"):
+            plm_files.append(file)
         else:
             otherfiles.append(file)
         
-    return confs, top_file, models, forcefield, index_file, data_file, dp_files, otherfiles
+    return confs, top_file, models, forcefield, index_file, data_file, dp_files, plm_files, otherfiles
 
 
 def log_ui():
-    logger.info('The task is displayed on "https://127.0.0.1:2746".')
-    logger.info('Artifacts (Files) are listed on "https://127.0.0.1:9001".')
+    if os.getenv("DFLOW_HOST") is not None:
+        logger.info('The task is displayed on %s.'%os.getenv("DFLOW_HOST"))
+    else:
+        logger.info('The task is displayed on "https://127.0.0.1:2746".')
+        logger.info('Artifacts (Files) are listed on "https://127.0.0.1:9001".')
 
 
 def main():
     args = parse_args()
     if args.command == "submit":
         logger.info("Preparing RiD ...")
-        confs, top_file, models, forcefield, index_file, data_file, dp_files, otherfiles = parse_submit(args)
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
         submit_rid(
             confs = confs,
             topology = top_file,
             rid_config = args.config,
             machine_config = args.machine,
+            workflow_id_defined = args.id,
             models = models,
             forcefield = forcefield,
             index_file = index_file,
@@ -299,13 +372,14 @@ def main():
         log_ui()
     elif args.command == "resubmit":
         logger.info("Preparing RiD ...")
-        confs, top_file, models, forcefield, index_file, data_file, dp_files, otherfiles = parse_submit(args)
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
         resubmit_rid(
             workflow_id=args.WORKFLOW_ID,
             confs = confs,
             topology = top_file,
             rid_config = args.config,
             machine_config = args.machine,
+            workflow_id_defined = args.id,
             iteration = args.iteration,
             pod = args.pod,
             models = models,
@@ -317,11 +391,23 @@ def main():
         )
         log_ui()
     elif args.command == "explore":
-        logger.info("RiD Exploration.")
-        return None
+        logger.info("RiD Exploration ..")
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
+        explore_rid(
+            confs = confs,
+            topology = top_file,
+            rid_config = args.config,
+            machine_config = args.machine,
+            models = models,
+            forcefield = forcefield,
+            index_file = index_file,
+            dp_files = dp_files,
+            otherfiles = otherfiles
+        )
+        log_ui()
     elif args.command == "label":
         logger.info("Labeling MD ...")
-        confs, top_file, models, forcefield, index_file, data_files, dp_files, otherfiles = parse_submit(args)
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
         label_rid(
             confs = confs,
             topology = top_file,
@@ -336,7 +422,7 @@ def main():
         log_ui()
     elif args.command == "relabel":
         logger.info("Labeling MD ...")
-        confs, top_file, models, forcefield, index_file, data_files, dp_files, otherfiles = parse_submit(args)
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
         relabel_rid(
             workflow_id=args.WORKFLOW_ID,
             confs = confs,
@@ -350,6 +436,48 @@ def main():
             otherfiles = otherfiles
         )
         log_ui()
+    elif args.command == "redim":
+        logger.info("Preparing MCMC ...")
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
+        redim_rid(
+            rid_config = args.config,
+            machine_config = args.machine,
+            models = models,
+            plm_out = plm_files,
+            workflow_id_defined = args.id,
+        )
+        log_ui()
+    elif args.command == "reredim":
+        logger.info("Preparing MCMC ...")
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
+        reredim_rid(
+            workflow_id=args.WORKFLOW_ID,
+            rid_config = args.config,
+            machine_config = args.machine,
+            workflow_id_defined = args.id,
+            models = models,
+            plm_out = plm_files,
+            pod = args.pod
+        )
+        log_ui()
+    elif args.command == "train":
+        logger.info("Training RiD FES ...")
+        confs, top_file, models, forcefield, index_file, data_file, dp_files,plm_files, otherfiles = parse_submit(args)
+        train_rid(
+            data = data_file,
+            rid_config = args.config,
+            machine_config = args.machine
+        )
+    elif args.command == "download":
+        logger.info("Downloading files ...")
+        download_rid(
+            workflow_id = args.WORKFLOW_ID,
+            pod = args.pod,
+            file = args.file,
+            iteration_start = args.iteration_start,
+            iteration_end = args.iteration_end,
+            outputs = args.outputs
+        )
     elif args.command == "port-forward":
         forward_ports()
     elif args.command == "ls":
diff --git a/rid/entrypoint/redim.py b/rid/entrypoint/redim.py
new file mode 100644
index 00000000..bc43299b
--- /dev/null
+++ b/rid/entrypoint/redim.py
@@ -0,0 +1,92 @@
+import json
+from pathlib import Path
+from typing import List, Union, Optional
+from rid.utils import load_json
+from copy import deepcopy
+import os
+
+from dflow import (
+    Workflow,
+    Step,
+    upload_artifact
+)
+
+from dflow.python import upload_packages
+from rid import SRC_ROOT
+upload_packages.append(SRC_ROOT)
+
+from rid.utils import normalize_resources
+from rid.superop.mcmc import MCMC
+from rid.op.mcmc_run import MCMCRun
+from rid.op.mcmc_plot import MCMCPlot
+
+def redim_rid(
+        rid_config: str,
+        machine_config: str,
+        models: Optional[Union[str, List[str]]] = None,
+        plm_out: Optional[Union[str, List[str]]] = None,
+        workflow_id_defined: Optional[str] = None
+    ):
+    with open(machine_config, "r") as mcg:
+        machine_config_dict = json.load(mcg)
+    resources = machine_config_dict["resources"]
+    tasks = machine_config_dict["tasks"]
+    normalized_resources = {}
+    for resource_type in resources.keys():
+        normalized_resources[resource_type] = normalize_resources(resources[resource_type])
+
+    mcmc_op = MCMC(
+        "mcmc",
+        MCMCRun,
+        MCMCPlot,
+        run_config = normalized_resources[tasks["mcmc_run_config"]],
+        plot_config = normalized_resources[tasks["mcmc_plot_config"]],
+        retry_times=1)
+    
+    jdata = deepcopy(load_json(rid_config))
+    mcmc_config = jdata["MCMC_Config"]
+    
+               
+    fe_models = []
+    assert isinstance(jdata["init_models"],list), "model input should be list."
+    for model in jdata["init_models"]:
+        fe_models.append(model)
+                
+    model_list = []
+            
+    if models is not None:
+        for model in models:
+            if os.path.basename(model) in fe_models:
+                model_list.append(model)
+
+    if len(model_list) == 0:
+        models_artifact = None
+    else:
+        models_artifact = upload_artifact([Path(p) for p in model_list], archive=None)
+        
+    plm_artifact =  None
+    if len(plm_out) != 0:
+        if isinstance(plm_out, str):
+            plm_artifact = upload_artifact(Path(plm_out), archive=None)
+        elif isinstance(plm_out, list):
+            plm_artifact = upload_artifact([Path(p) for p in plm_out], archive=None)
+    
+    task_names = []
+    for index in range(len(model_list)):
+        task_names.append("%03d"%index)
+
+    rid_steps = Step("rid-mcmc",
+            mcmc_op,
+            artifacts={
+                "models": models_artifact,
+                "plm_out": plm_artifact
+            },
+            parameters={
+                "mcmc_config": mcmc_config,
+                "task_names": task_names,
+                "block_tag" : "000"
+            },
+        )
+    wf = Workflow("rid-mcmc", pod_gc_strategy="OnPodSuccess", parallelism=10, id = workflow_id_defined)
+    wf.add(rid_steps)
+    wf.submit()
diff --git a/rid/entrypoint/reredim.py b/rid/entrypoint/reredim.py
new file mode 100644
index 00000000..035e3905
--- /dev/null
+++ b/rid/entrypoint/reredim.py
@@ -0,0 +1,120 @@
+import json
+from pathlib import Path
+from typing import List, Union, Optional
+from rid.utils import load_json
+from copy import deepcopy
+import os
+
+from dflow import (
+    Workflow,
+    Step,
+    upload_artifact
+)
+
+from dflow.python import upload_packages
+from rid import SRC_ROOT
+upload_packages.append(SRC_ROOT)
+
+from rid.utils import normalize_resources
+from rid.superop.mcmc import MCMC
+from rid.op.mcmc_run import MCMCRun
+from rid.op.mcmc_plot import MCMCPlot
+
+def reredim_rid(
+        workflow_id: str,
+        rid_config: str,
+        machine_config: str,
+        models: Optional[Union[str, List[str]]] = None,
+        plm_out: Optional[Union[str, List[str]]] = None,
+        workflow_id_defined: Optional[str] = None,
+        pod: Optional[str] = None
+    ):
+    with open(machine_config, "r") as mcg:
+        machine_config_dict = json.load(mcg)
+    resources = machine_config_dict["resources"]
+    tasks = machine_config_dict["tasks"]
+    normalized_resources = {}
+    for resource_type in resources.keys():
+        normalized_resources[resource_type] = normalize_resources(resources[resource_type])
+
+    mcmc_op = MCMC(
+        "mcmc",
+        MCMCRun,
+        MCMCPlot,
+        run_config = normalized_resources[tasks["mcmc_run_config"]],
+        plot_config = normalized_resources[tasks["mcmc_plot_config"]],
+        retry_times=None)
+    
+    jdata = deepcopy(load_json(rid_config))
+    mcmc_config = jdata["MCMC_Config"]
+    
+               
+    fe_models = []
+    assert isinstance(jdata["init_models"],list), "model input should be list."
+    for model in jdata["init_models"]:
+        fe_models.append(model)
+                
+    model_list = []
+            
+    if models is not None:
+        for model in models:
+            if os.path.basename(model) in fe_models:
+                model_list.append(model)
+
+    if len(model_list) == 0:
+        models_artifact = None
+    else:
+        models_artifact = upload_artifact([Path(p) for p in model_list], archive=None)
+        
+    plm_artifact =  None
+    if len(plm_out) != 0:
+        if isinstance(plm_out, str):
+            plm_artifact = upload_artifact(Path(plm_out), archive=None)
+        elif isinstance(plm_out, list):
+            plm_artifact = upload_artifact([Path(p) for p in plm_out], archive=None)
+    
+    
+    task_names = []
+    for index in range(len(model_list)):
+        task_names.append("%03d"%index)
+
+    rid_steps = Step("rid-mcmc",
+            mcmc_op,
+            artifacts={
+                "models": models_artifact,
+                "plm_out": plm_artifact
+            },
+            parameters={
+                "mcmc_config": mcmc_config,
+                "task_names": task_names,
+                "block_tag" : "000"
+            },
+        )
+    
+    old_workflow = Workflow(id=workflow_id)
+    all_steps = old_workflow.query_step()
+
+    succeeded_steps = []
+    restart_flag = 1
+    for step in all_steps:
+        if step["type"] == "Pod":
+            pod_key = step["key"]
+            if pod_key is not None:
+                pod_key_list = pod_key.split("-")
+                pod_step_1 = "-".join(pod_key_list[1:-1])
+                pod_step_2 = "-".join(pod_key_list[1:])
+                if pod is not None:
+                    if pod_step_1 == pod or pod_step_2 == pod:
+                        restart_flag = 0
+                else:
+                    if step["phase"] != "Succeeded":
+                        restart_flag = 0
+                    else:
+                        restart_flag = 1
+            
+            if restart_flag == 1:
+                succeeded_steps.append(step)
+                
+    wf = Workflow("rid-mcmc-continue", pod_gc_strategy="OnPodSuccess", parallelism=10, id = workflow_id_defined)
+    wf.add(rid_steps)
+    wf.submit(reuse_step=succeeded_steps)
diff --git a/rid/entrypoint/resubmit.py b/rid/entrypoint/resubmit.py
index f6fe7944..a1df973e 100644
--- a/rid/entrypoint/resubmit.py
+++ b/rid/entrypoint/resubmit.py
@@ -24,6 +24,7 @@ def resubmit_rid(
         topology: Optional[str],
         rid_config: str,
         machine_config: str,
+        workflow_id_defined: Optional[str] = None,
         iteration: Optional[str] = None,
         pod: Optional[str] = None,
         models: Optional[Union[str, List[str]]] = None,
@@ -50,6 +51,7 @@ def resubmit_rid(
         run_select_config = normalized_resources[tasks["run_select_config"]],
         prep_data_config = normalized_resources[tasks["prep_data_config"]],
         run_train_config = normalized_resources[tasks["run_train_config"]],
+        model_devi_config = normalized_resources[tasks["model_devi_config"]],
         workflow_steps_config = normalized_resources[tasks["workflow_steps_config"]],
         retry_times=None
     )
@@ -181,10 +183,11 @@ def resubmit_rid(
                 if pod_key is not None:
                     pod_key_list = pod_key.split("-")
                     pod_iter = int(pod_key_list[1])
-                    pod_step = "-".join(pod_key_list[2:-1])
+                    pod_step_1 = "-".join(pod_key_list[2:-1])
+                    pod_step_2 = "-".join(pod_key_list[2:])
                     if iteration is not None:
                         if pod is not None:
-                            if pod_iter == int(iteration) and pod_step == pod:
+                            if pod_iter == int(iteration) and (pod_step_1 == pod or pod_step_2 == pod):
                                 restart_flag = 0
                         else:
                             if pod_iter == int(iteration):
@@ -197,6 +200,6 @@ def resubmit_rid(
                 
                 if restart_flag == 1:
                     succeeded_steps.append(step)
-    wf = Workflow("reinforced-dynamics-continue", pod_gc_strategy="OnPodSuccess", parallelism=50)
+    wf = Workflow("reinforced-dynamics-continue", pod_gc_strategy="OnPodSuccess", parallelism=50, id = workflow_id_defined)
     wf.add(rid_steps)
     wf.submit(reuse_step=succeeded_steps)
diff --git a/rid/entrypoint/server.py b/rid/entrypoint/server.py
index 666e9e85..781a61d7 100644
--- a/rid/entrypoint/server.py
+++ b/rid/entrypoint/server.py
@@ -41,16 +41,16 @@ def forward_ports():
         logger.info('Port "agro-server" has been launched and running.')
     else:
         port = PORTS["argo"]
-        os.system(f"nohup kubectl -n {argo_namespace} port-forward deployment/argo-server {port}:{port} --address 0.0.0.0 > /tmp/argo-server.out 2>&1 &")
+        os.system(f"nohup minikube kubectl -- -n {argo_namespace} port-forward deployment/argo-server {port}:{port} --address 0.0.0.0 > /tmp/argo-server.out 2>&1 &")
     
     if port_status["minio-server"] == 1:
         logger.info('Port "minio-server" has been launched and running.')
     else:
         port = PORTS["minio-server"]
-        os.system(f"nohup kubectl -n {argo_namespace} port-forward deployment/minio {port}:{port} --address 0.0.0.0 > /tmp/minio-server.out 2>&1 &")
+        os.system(f"nohup minikube kubectl -- -n {argo_namespace} port-forward deployment/minio {port}:{port} --address 0.0.0.0 > /tmp/minio-server.out 2>&1 &")
     
     if port_status["minio-ui"] == 1:
         logger.info('Port "minio-ui" has been launched and running.')
     else:
         port = PORTS["minio-ui"]
-        os.system(f"nohup kubectl -n {argo_namespace} port-forward deployment/minio {port}:{port} --address 0.0.0.0 > /tmp/minio-ui.out 2>&1 &")
+        os.system(f"nohup minikube kubectl -- -n {argo_namespace} port-forward deployment/minio {port}:{port} --address 0.0.0.0 > /tmp/minio-ui.out 2>&1 &")
diff --git a/rid/entrypoint/submit.py b/rid/entrypoint/submit.py
index 618993e1..ab1b0c82 100644
--- a/rid/entrypoint/submit.py
+++ b/rid/entrypoint/submit.py
@@ -22,6 +22,7 @@
 from rid.superop.label import Label
 from rid.op.prep_label import PrepLabel, CheckLabelInputs
 from rid.op.run_label import RunLabel
+from rid.op.label_stats import LabelStats
 from rid.superop.selector import Selector
 from rid.op.prep_select import PrepSelect
 from rid.op.run_select import RunSelect
@@ -29,6 +30,7 @@
 from rid.op.prep_data import CollectData, MergeData
 from rid.superop.blocks import IterBlock, InitBlock
 from rid.op.run_train import TrainModel
+from rid.op.run_model_devi import RunModelDevi
 from rid.op.adjust_trust_level import AdjustTrustLevel
 from rid.flow.loop import ReinforcedDynamics
 
@@ -42,6 +44,7 @@ def prep_rid_op(
     run_select_config,
     prep_data_config,
     run_train_config,
+    model_devi_config,
     workflow_steps_config,
     retry_times
     ):
@@ -59,6 +62,7 @@ def prep_rid_op(
         CheckLabelInputs,
         PrepLabel,
         RunLabel,
+        LabelStats,
         prep_label_config,
         run_label_config,
         retry_times=retry_times)
@@ -85,7 +89,9 @@ def prep_rid_op(
         label_op,
         data_op,
         TrainModel,
+        RunModelDevi,
         run_train_config,
+        model_devi_config,
         retry_times=retry_times
     )
 
@@ -97,8 +103,10 @@ def prep_rid_op(
         data_op,
         AdjustTrustLevel,
         TrainModel,
+        RunModelDevi,
         workflow_steps_config,
         run_train_config,
+        model_devi_config,
         retry_times=retry_times)
     
     rid_op = ReinforcedDynamics(
@@ -114,6 +122,7 @@ def submit_rid(
         topology: Optional[str],
         rid_config: str,
         machine_config: str,
+        workflow_id_defined: Optional[str] = None,
         models: Optional[Union[str, List[str]]] = None,
         forcefield: Optional[str] = None,
         index_file: Optional[str] = None,
@@ -138,6 +147,7 @@ def submit_rid(
         run_select_config = normalized_resources[tasks["run_select_config"]],
         prep_data_config = normalized_resources[tasks["prep_data_config"]],
         run_train_config = normalized_resources[tasks["run_train_config"]],
+        model_devi_config = normalized_resources[tasks["model_devi_config"]],
         workflow_steps_config = normalized_resources[tasks["workflow_steps_config"]],
         retry_times=1
     )
@@ -255,6 +265,6 @@ def submit_rid(
             },
             parameters={}
         )
-    wf = Workflow("reinforced-dynamics", pod_gc_strategy="OnPodSuccess", parallelism=50)
+    wf = Workflow("reinforced-dynamics", pod_gc_strategy="OnPodSuccess", parallelism=50, id = workflow_id_defined)
     wf.add(rid_steps)
     wf.submit()
diff --git a/rid/entrypoint/train.py b/rid/entrypoint/train.py
new file mode 100644
index 00000000..9c6ece94
--- /dev/null
+++ b/rid/entrypoint/train.py
@@ -0,0 +1,86 @@
+import json
+from pathlib import Path
+from rid.utils import load_json
+from copy import deepcopy
+
+from dflow import (
+    Workflow,
+    Step,
+    Steps,
+    upload_artifact,
+    argo_range,
+    argo_len
+)
+
+from dflow.python import(
+    PythonOPTemplate,
+    upload_packages,
+    Slices,
+)
+
+from rid import SRC_ROOT
+upload_packages.append(SRC_ROOT)
+
+from rid.utils import normalize_resources,init_executor
+from rid.op.run_train import TrainModel
+from rid.constants import model_tag_fmt
+
+
+def train_rid(
+        data: str,
+        rid_config: str,
+        machine_config: str
+    ):
+    with open(machine_config, "r") as mcg:
+        machine_config_dict = json.load(mcg)
+    resources = machine_config_dict["resources"]
+    tasks = machine_config_dict["tasks"]
+    normalized_resources = {}
+    for resource_type in resources.keys():
+        normalized_resources[resource_type] = normalize_resources(resources[resource_type])
+
+    jdata = deepcopy(load_json(rid_config))
+    cv_config = jdata["CV"]
+    train_config = jdata["Train"]
+    numb_models = train_config["numb_models"]
+    angular_mask = cv_config["angular_mask"]
+    
+    run_train_config = normalized_resources[tasks["run_train_config"]]
+    run_train_config = deepcopy(run_train_config)
+    train_template_config = run_train_config.pop('template_config')
+    train_executor = init_executor(run_train_config.pop('executor'))
+    
+    model_tags = []
+    for idx in range(numb_models):
+        model_tags.append(model_tag_fmt.format(idx=idx))
+        
+    data_artifact = upload_artifact(Path(data), archive=None)
+    
+    train_op = TrainModel
+    train_steps=Step(
+        "train",
+        template=PythonOPTemplate(
+            train_op,
+            python_packages = None,
+            retry_on_transient_error = 1,
+            slices=Slices("{{item}}",
+                input_parameter=["model_tag"],
+                output_artifact=["model","train_fig"]),
+            **train_template_config,
+        ),
+        parameters={
+            "model_tag": model_tags,
+            "angular_mask": angular_mask,
+            "train_config": train_config,
+        },
+        artifacts={
+            "data": data_artifact
+        },
+        executor = train_executor,
+        with_param = argo_range(len(model_tags)),
+        key = "run-train"+"-{{item}}",
+        **run_train_config,
+    )
+    wf = Workflow("rid-train", pod_gc_strategy="OnPodSuccess", parallelism=10)
+    wf.add(train_steps)
+    wf.submit()
diff --git a/rid/flow/loop.py b/rid/flow/loop.py
index 52309abf..008f676d 100644
--- a/rid/flow/loop.py
+++ b/rid/flow/loop.py
@@ -45,12 +45,14 @@ def __init__(
             "weights": InputParameter(type=Optional[Union[np.ndarray, List]]),
             "max_selection": InputParameter(type=int),
             "numb_cluster_threshold": InputParameter(type=float, value=30),
+            "std_threshold": InputParameter(type=float, value=5.0),
             "dt": InputParameter(type=float, value=0.02),
             "output_freq": InputParameter(type=float, value=2500),
             "slice_mode": InputParameter(type=str, value="gmx"),
             "label_config": InputParameter(type=Dict),
             "tail": InputParameter(type=float, value=0.9),
             "train_config": InputParameter(type=Dict),
+            "type_map": InputParameter(type=List, value = []),
             "adjust_amplifier": InputParameter(type=float, value=1.5),
             "max_level_multiple": InputParameter(type=float, value=8.0),
         }
@@ -168,9 +170,11 @@ def _loop (
             "weights": steps.inputs.parameters["weights"],
             "max_selection": steps.inputs.parameters["max_selection"],
             "numb_cluster_threshold": steps.inputs.parameters["numb_cluster_threshold"],
+            "std_threshold": steps.inputs.parameters["std_threshold"],
             "dt": steps.inputs.parameters["dt"],
             "output_freq": steps.inputs.parameters["output_freq"],
             "slice_mode": steps.inputs.parameters["slice_mode"],
+            "type_map": steps.inputs.parameters["type_map"],
             "label_config": steps.inputs.parameters["label_config"],
             "train_config": steps.inputs.parameters["train_config"]
         },
@@ -209,9 +213,11 @@ def _loop (
             "weights": steps.inputs.parameters["weights"],
             "max_selection": steps.inputs.parameters["max_selection"],
             "numb_cluster_threshold": steps.inputs.parameters["numb_cluster_threshold"],
+            "std_threshold": steps.inputs.parameters["std_threshold"],
             "dt": steps.inputs.parameters["dt"],
             "output_freq": steps.inputs.parameters["output_freq"],
             "slice_mode": steps.inputs.parameters["slice_mode"],
+            "type_map": steps.inputs.parameters["type_map"],
             "label_config": steps.inputs.parameters["label_config"],
             "train_config": steps.inputs.parameters["train_config"]
         },
@@ -408,9 +414,11 @@ def _rid(
             "angular_mask": prep_rid.outputs.parameters["angular_mask"],
             "weights": prep_rid.outputs.parameters["weights"],
             "max_selection": prep_rid.outputs.parameters["max_selection"],
+            "std_threshold": prep_rid.outputs.parameters["std_threshold"],
             "dt": prep_rid.outputs.parameters["dt"],
             "output_freq": prep_rid.outputs.parameters["output_freq"],
             "slice_mode": prep_rid.outputs.parameters["slice_mode"],
+            "type_map": prep_rid.outputs.parameters["type_map"],
             "label_config": prep_rid.outputs.parameters["label_config"],
             "train_config": prep_rid.outputs.parameters["train_config"]
         },
@@ -454,9 +462,11 @@ def _rid(
             "weights": prep_rid.outputs.parameters["weights"],
             "max_selection": prep_rid.outputs.parameters["max_selection"],
             "numb_cluster_threshold": prep_rid.outputs.parameters["numb_cluster_threshold"],
+            "std_threshold": prep_rid.outputs.parameters["std_threshold"],
             "dt": prep_rid.outputs.parameters["dt"],
             "output_freq": prep_rid.outputs.parameters["output_freq"],
             "slice_mode": prep_rid.outputs.parameters["slice_mode"],
+            "type_map": prep_rid.outputs.parameters["type_map"],
             "label_config": prep_rid.outputs.parameters["label_config"],
             "train_config": prep_rid.outputs.parameters["train_config"]
         },
diff --git a/rid/mcmc/__init__.py b/rid/mcmc/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/rid/mcmc/walker.py b/rid/mcmc/walker.py
new file mode 100644
index 00000000..1efcc504
--- /dev/null
+++ b/rid/mcmc/walker.py
@@ -0,0 +1,140 @@
+import numpy as np
+import tensorflow as tf
+
+# kinetic enery in eV
+kbT = (8.617343E-5) * 300 
+beta = 1.0 / kbT
+
+# ev to kj/mol
+f_cvt = 96.485
+
+def distance(cv1,cv2):
+    d = np.sum((cv1 - cv2)**2,axis=1)
+    return d
+
+def s(cv, lm, xlist):
+    s = 0
+    n = 1e-10
+    for i in range(1,10):
+       s +=  i*np.exp(-lm*distance(cv,xlist[i-1]))
+       n += np.exp(-lm*distance(cv,xlist[i-1]))
+    return s/n
+
+def z(cv, lm, xlist):
+    z = 0
+    for i in range(9):
+        z += np.exp(-lm*distance(cv,xlist[i]))
+    p = np.ones(shape=z.shape)*1e-10
+    z = np.log(np.max((z,p),axis=0))
+    return -z/lm
+
+class Walker(object):
+    def __init__(self, fd, nw, sess, type, cv_lower, cv_upper):
+        self.type = type
+        self._full_dim = fd
+        self._num_walker = nw
+        self._move_scale = 0.5
+        self._sess = sess
+        self._shape = (self._num_walker, self._full_dim)
+        self._cv_range = cv_upper
+        self._cv_shift = cv_lower
+        # absolute coordinate
+        self._position = np.random.uniform(size=self._shape)*(np.array(self._cv_range) - np.array(self._cv_shift)) + np.array(self._cv_shift)
+        self._energy = np.zeros([self._num_walker])
+        self._force = np.zeros([self._num_walker, self._full_dim])
+        self._sample_step = 20
+        self._acp_ratio_lb = 0.15
+        self._acp_ratio_ub = 0.75
+        self.max_scale = 5
+        self.min_scale = 0.01
+        self.inc_scale_fac = 1.25
+
+    def sample(self, compute_ef, inter_step=1):
+        acp_ratio = []
+        self._energy, self._force = compute_ef(self._sess, self._position)
+        for _ in range(inter_step):
+            if self.type == "dih":
+                position_new = np.mod(self._position + np.random.normal(scale=self._move_scale,size=self._shape), 2*np.pi)
+            elif self.type == "dis":
+                position_new = self._position + np.random.normal(scale=self._move_scale*0.08, size=self._shape)
+                indices = np.where(position_new < 0.1)
+                position_new[indices]=0.13 + np.random.normal(scale=0.01, size=indices[0].shape)
+                indices = np.where(position_new > 9.9)
+                position_new[indices]=9.87 + np.random.normal(scale=0.01, size=indices[0].shape)
+            else:
+                raise ValueError("Undefined cv type, only support 'dih' and 'dis' type")
+            energy_new, force_new = compute_ef(self._sess, position_new)
+            # in case of overflow
+            prob_ratio = np.exp(np.minimum(- beta * (energy_new - self._energy), 0))
+            idx = np.random.uniform(size=self._num_walker) < np.reshape(prob_ratio, [-1])
+
+            self._position[idx, :] = position_new[idx, :]
+            self._energy[idx] = energy_new[idx]
+            self._force[idx] = force_new[idx]
+            acp_ratio.append(np.mean(idx))
+            # elif self.type == "dis":
+            #     for ii in range(self._num_walker):
+            #         if (position_new[ii,:] > self._cv_shift).all() and (position_new[ii,:] < self._cv_range).all() and idx[ii]:
+            #             self._position[ii, :] = position_new[ii, :]
+            #             self._energy[ii] = energy_new[ii]
+            #             self._force[ii] = force_new[ii]
+            #             acp_ratio.append(np.mean(idx))
+            # else:
+            #     raise ValueError("Undefined cv type, only support 'dih' and 'dis' type")
+        acp_ratio = np.mean(acp_ratio)
+        if acp_ratio > self._acp_ratio_ub:
+            # move_scale is too small
+            self._move_scale = np.min((self._move_scale*self.inc_scale_fac,self.max_scale),axis=0)
+            print(
+                "Increase move_scale to %s due to high acceptance ratio: %f" % (
+                    self._move_scale, acp_ratio))
+            # print(self._position[:5, :, :])
+        elif acp_ratio < self._acp_ratio_lb:
+            # move_scale is too large
+            self._move_scale = np.max((self._move_scale/self.inc_scale_fac,self.min_scale),axis=0)
+            print(
+                "Decrease move_scale to %s due to low acceptance ratio: %f" % (
+                    self._move_scale, acp_ratio))
+        return self._position, self._energy, self._force
+
+# project on all 1d points
+def my_hist1d(pp, xx, delta, fd):
+    my_hist = np.zeros((fd, len(xx)))
+    for ii in range(pp.shape[0]):   ###trj_num
+        for jj in range(fd):        ###cv_num
+            my_hist[jj, np.int(pp[ii,jj]//delta)] += 1
+    my_hist /= (pp.shape[0] * delta)
+    return my_hist
+
+# project on specific 2d point
+def my_hist2d(pp, xx, yy, delta, cv1, cv2):
+    my_hist = np.zeros((1, len(xx), len(yy)))
+    for ii in range(pp.shape[0]):
+        my_hist[0, np.int(pp[ii,cv1]//delta), np.int(pp[ii,cv2]//delta)] += 1
+    my_hist /= (pp.shape[0] * delta * delta)
+    return my_hist
+
+# project on all path CV 1d
+def my_hist1d_path(pp, xx, delta, lm, xlist, proj_index):
+    my_hist = np.zeros((2, len(xx)))
+    cv = pp[:,proj_index]
+    slist = s(cv, lm, xlist)
+    zlist = z(cv, lm, xlist)
+    for ii in range(cv.shape[0]):
+        if 0 < slist[ii] < 10 and -10 < zlist[ii] < 0:
+            my_hist[0, np.int((slist[ii])//delta)] += 1
+            my_hist[1, np.int((zlist[ii]+10)//delta)] += 1
+    my_hist /= (pp.shape[0] * delta)
+    return my_hist
+
+# project on path CV 2d
+def my_hist2d_path(pp, xx, yy, delta, lm, xlist, proj_index):
+    my_hist = np.zeros((1, len(xx), len(yy)))
+    cv = pp[:,proj_index]
+    slist = s(cv, lm, xlist)
+    zlist = z(cv, lm, xlist)
+    for ii in range(cv.shape[0]):
+        if 0 < slist[ii] < 10 and -10 < zlist[ii] < 0:
+            my_hist[0, np.int((slist[ii])//delta), np.int((zlist[ii]+10)//delta)] += 1
+    my_hist /= (pp.shape[0] * delta * delta)
+    return my_hist
\ No newline at end of file
diff --git a/rid/op/label_stats.py b/rid/op/label_stats.py
new file mode 100644
index 00000000..53b7e3a7
--- /dev/null
+++ b/rid/op/label_stats.py
@@ -0,0 +1,116 @@
+from typing import Dict, List
+from pathlib import Path
+from dflow.python import (
+    OP,
+    OPIO,
+    OPIOSign,
+    Artifact,
+    Parameter,
+    BigParameter
+)
+import numpy as np
+from rid.utils import set_directory
+from matplotlib import pyplot as plt
+from rid.constants import mf_std_fig
+import os
+
+
+class LabelStats(OP):
+
+    r"""LabelStats OP calculate the std of all the labeling steps, remove the ones with high std.
+    """
+
+    @classmethod
+    def get_input_sign(cls):
+        return OPIOSign(
+            {
+                "cv_forces": Artifact(List[Path]),
+                "mf_info": Artifact(List[Path]), 
+                "std_threshold": float  
+            }
+        )
+
+    @classmethod
+    def get_output_sign(cls):
+        return OPIOSign(
+            {
+                "mf_std_fig":  Artifact(Path, archive=None),
+                "cv_forces": Artifact(List[Path], archive=None)
+            }
+        )
+
+    @OP.exec_sign_check
+    def execute(
+        self,
+        op_in: OPIO,
+    ) -> OPIO:
+        r"""Execute the OP.
+
+        Parameters
+        ----------
+        op_in : dict
+            Input dict with components:
+        
+            - "cv_forces": (`Artifact(List[Path])`) 
+            - "mf_info": (`Artifact(List[Path])`) 
+            - "std_threshold": (`Float`)
+
+        Returns
+        -------
+            Output dict with components:
+        
+            - "mf_std_fig":  (`Artifact(Path)`)
+            - "cv_forces":(`Artifact(List[Path])`) 
+            
+        """
+        assert len(op_in["mf_info"]) == len(op_in["cv_forces"])
+        cv_forces_list = [_ for _ in op_in["cv_forces"] if _ is not None]
+        cv_forces_list = np.array(cv_forces_list)
+        cv_force_file = cv_forces_list[0]
+        cv_force = np.loadtxt(cv_force_file)
+        cv_dim = int(cv_force.shape[0]//2)
+        # extract the std from the mf_info
+        mf_all_std_list = []
+        for mf_info in op_in["mf_info"]:
+            if mf_info and os.path.exists(mf_info):
+                with open(mf_info) as f:
+                    while True:
+                        line = f.readline()
+                        line_list = line.strip().split(" ")[:3]
+                        if line_list == ['mean', 'force', 'std']:
+                            break
+                    mf_std_line = line.strip().split(" ")[7:]
+                    mf_std_list = [float(i) for i in mf_std_line]
+                    mf_all_std_list.append(mf_std_list)
+        mf_all_std_list = np.array(mf_all_std_list)
+        
+        higher_index = set()
+        for i, row in enumerate(mf_all_std_list):
+            for j, num in enumerate(row):
+                if num > op_in["std_threshold"]:                    
+                    higher_index.add(i)
+        higher_index_list = list(higher_index)
+        print("higher index list", list(cv_forces_list[higher_index_list]))
+        mf_all_std_list_modified = np.delete(mf_all_std_list, higher_index_list, axis=0)
+        cv_forces_list_modified = np.delete(cv_forces_list, higher_index_list, axis=0)
+        assert len(mf_all_std_list_modified) == len(cv_forces_list_modified)
+        
+        task_path = Path("label_std")
+        task_path.mkdir(exist_ok=True, parents=True)
+        with set_directory(task_path):
+            plt.figure(figsize=(8,6))
+            for cv_index in range(cv_dim):
+                plt.hist(mf_all_std_list_modified[:,cv_index], 100, label = "cv%s"%(cv_index+1))
+            plt.legend()
+            plt.xlabel("mean force std")
+            plt.ylabel("frequency")
+            plt.savefig(mf_std_fig)
+             
+            
+        op_out = OPIO({
+            "mf_std_fig":  task_path.joinpath(mf_std_fig),
+            "cv_forces": list(cv_forces_list_modified)
+        })
+        return op_out
+
+    
\ No newline at end of file
diff --git a/rid/op/mcmc_plot.py b/rid/op/mcmc_plot.py
new file mode 100644
index 00000000..8f17bdc7
--- /dev/null
+++ b/rid/op/mcmc_plot.py
@@ -0,0 +1,259 @@
+from typing import List, Optional, Dict
+from pathlib import Path
+import numpy as np
+from dflow.python import (
+    OP,
+    OPIO,
+    OPIOSign,
+    Artifact,
+    Parameter,
+    BigParameter
+)
+from rid.utils import save_txt, set_directory
+from rid.constants import mcmc_1cv_dir_name, mcmc_1cv_name, mcmc_2cv_name,mcmc_1cv_fig, mcmc_2cv_fig, mcmc_2cv_fig_separate
+from matplotlib import pyplot as plt
+import os
+
+class MCMCPlot(OP):
+    """
+    `MCMC_Plot` plot the reduced free energy surface produced by MCMC run.
+    """
+
+    @classmethod
+    def get_input_sign(cls):
+        return OPIOSign(
+            {
+                "mcmc_1cv": Artifact(List[Path]),
+                "mcmc_2cv": Artifact(List[Path]),
+                "plm_out": Artifact(List[Path], optional=True),
+                "mcmc_config": BigParameter(Dict)
+            }
+        )
+
+    @classmethod
+    def get_output_sign(cls):
+        return OPIOSign(
+            {
+                "mcmc_fig": Artifact(Path, archive = None)
+            }
+        )
+
+    @OP.exec_sign_check
+    def execute(
+        self,
+        op_in: OPIO,
+    ) -> OPIO:
+
+        r"""Execute the OP.
+        
+        Parameters
+        ----------
+        op_in : dict
+            Input dict with components:
+
+            - `mcmc_config`: Dict,
+            - `mcmc_cv`: Artifact(List[Path]),
+          
+        Returns
+        -------
+            Output dict with components:
+        
+            - `mcmc_fig`: (`Artifact(Path)`)
+        """
+        mcmc_config = op_in["mcmc_config"]
+        cv_type = mcmc_config["cv_type"]
+        bins = mcmc_config["bins"]
+        cv_dim = mcmc_config["cv_dimension"]
+        proj_info = mcmc_config["proj_info"]
+        proj_mode = proj_info["proj_mode"]
+        
+        if proj_mode == "cv":
+            proj_1d_num = cv_dim
+        elif proj_mode == "path":
+            proj_1d_num = 2
+        else:
+            raise ValueError("Invalid cv type")
+        
+        task_path = Path("mcmc_fig")
+        task_path.mkdir(exist_ok=True, parents=True)
+        with set_directory(task_path):
+            fourdata_1cv = []
+            for dir in op_in["mcmc_1cv"]:
+                all_1cv = []
+                for cv_index in range(proj_1d_num):
+                    file = dir/mcmc_1cv_name.format(tag = cv_index)
+                    pmf = np.loadtxt(file)
+                    pmf = pmf - np.min(pmf)
+                    allda=[]
+                    proj_iterations = pmf.shape[0]
+                    if cv_type == "dih":
+                        xedges = np.linspace(0, 2*np.pi, bins)
+                        for i in range(bins):
+                            temp = []
+                            # this is to change to dihedral dimension to (-pi,pi) which corresponds to the gmx output
+                            if xedges[i]>=np.pi:
+                                temp.append(xedges[i] - np.pi*2)
+                            else:
+                                temp.append(xedges[i] + np.pi*2/(bins-1))
+                            if i == (bins - 1):
+                                temp.append(pmf[:,0])
+                            else:
+                                temp.append(pmf[:,i])
+                            temp = np.array(temp)
+                            allda.append(temp)
+                    elif cv_type == "dis":
+                        xedges = np.linspace(0, 10, bins)
+                        for i in range(bins):
+                            temp = []
+                            temp.append(xedges[i])
+                            temp.append(pmf[:,i])
+                            temp = np.array(temp)
+                            allda.append(temp)
+                    newarray=np.array(allda)
+                    idex=np.argsort(newarray[:,0])
+                    sorted_data = newarray[idex,:]
+                    print("sorted data shape",sorted_data.shape)
+                    sorted_data = np.stack(sorted_data[:,1])
+                    all_1cv.append(sorted_data)
+                all_1cv = np.array(all_1cv)
+                print("all_1cv shape", all_1cv.shape)
+                fourdata_1cv.append(all_1cv)
+            avedata_1cv = np.mean(np.array(fourdata_1cv),axis=0)
+            print("avedata shape", avedata_1cv.shape)
+            # make 1cv plot
+            if not os.path.exists(mcmc_1cv_dir_name):
+                os.makedirs(mcmc_1cv_dir_name)
+
+            for cv_index in range(proj_1d_num):
+                if proj_mode == "cv":
+                    if cv_type == "dih":
+                        xedges = np.linspace(-np.pi, np.pi, bins)
+                    elif cv_type == "dis":
+                        xedges = np.linspace(0, 10, bins)
+                    plt.figure(figsize=(8, 6))
+                    for proj_iter in range(proj_iterations):
+                        plt.plot(xedges,avedata_1cv[cv_index][:,proj_iter], label = proj_iter)
+                    plt.xlabel(r'cv')
+                    plt.ylabel(r'free energy (kcal/mol)')
+                    plt.legend()
+                    plt.savefig(mcmc_1cv_dir_name+"/"+mcmc_1cv_fig.format(tag = cv_index),dpi=600,bbox_inches='tight')
+                elif proj_mode == "path":
+                    path_edges = [np.linspace(0, 10, bins),np.linspace(-10, 0, bins)]
+                    plt.figure(figsize=(8, 6))
+                    for proj_iter in range(proj_iterations):
+                        plt.plot(path_edges[cv_index],avedata_1cv[cv_index][:,proj_iter], label = proj_iter)
+                    plt.xlabel(r'cv')
+                    plt.ylabel(r'free energy (kcal/mol)')
+                    plt.legend()
+                    plt.savefig(mcmc_1cv_dir_name+"/"+mcmc_1cv_fig.format(tag = cv_index),dpi=600,bbox_inches='tight')
+                    
+            fourdata_2cv = []
+            for file in op_in["mcmc_2cv"]:
+                pmf = np.loadtxt(file)
+                pmf = pmf - np.min(pmf)
+                allda=[]
+                if cv_type == "dih":
+                    xedges = np.linspace(0, 2*np.pi, bins)
+                    yedges = np.linspace(0, 2*np.pi, bins)
+                    for i in range(bins):
+                        for j in range(bins):
+                            temp = []
+                            # this is to change to dihedral dimension to (-pi,pi) which corresponds to the gmx output
+                            if xedges[i]>=np.pi:
+                                temp.append(xedges[i] - np.pi*2)
+                            else:
+                                temp.append(xedges[i] + np.pi*2/(bins-1))
+                            if yedges[j]>=np.pi:
+                                temp.append(yedges[j] - np.pi*2)
+                            else:
+                                temp.append(yedges[j] + np.pi*2/(bins-1))
+                            if j == (bins -1):
+                                temp.append(pmf[i][0])
+                            elif i == (bins - 1):
+                                temp.append(pmf[0][j])
+                            else:
+                                temp.append(pmf[i][j])
+                            temp = np.array(temp)
+                            allda.append(temp)
+                elif cv_type == "dis":
+                    xedges = np.linspace(0, 10, bins)
+                    yedges = np.linspace(0, 10, bins)
+                    for i in range(bins):
+                        for j in range(bins):
+                            temp = []
+                            temp.append(xedges[i])
+                            temp.append(yedges[j])
+                            temp.append(pmf[i][j])
+                            temp = np.array(temp)
+                            allda.append(temp)
+                newarray=np.array(allda)
+                idex=np.lexsort([newarray[:,1], newarray[:,0]])
+                sorted_data = newarray[idex,:]
+                fourdata_2cv.append(sorted_data[:,2])
+            avedata_2cv = np.mean(np.array(fourdata_2cv),axis=0)    
+            # make 2cv plot
+            if cv_type == "dih":
+                xedges = np.linspace(-np.pi, np.pi, bins)
+                yedges = np.linspace(-np.pi, np.pi, bins)
+            elif cv_type == "dis":
+                if proj_mode == "cv":
+                    xedges = np.linspace(0, 10, bins)
+                    yedges = np.linspace(0, 10, bins)
+                elif proj_mode == "path":
+                    xedges = np.linspace(0, 10, bins)
+                    yedges = np.linspace(-10, 0, bins)
+            plt.figure(figsize=(8, 6))
+            cmap = plt.cm.get_cmap("jet_r")
+            # Define percentiles for the levels
+            upper_perc = np.percentile(np.unique(fourdata_2cv[0]), 99)
+            CS = plt.contourf(xedges,yedges,avedata_2cv.reshape(bins,bins),levels = np.linspace(0,upper_perc,101),cmap=cmap,extend="max")
+            
+            if op_in["plm_out"] is not None:
+                for cv_output in op_in["plm_out"]:
+                    cv_point = np.loadtxt(cv_output)
+                    assert cv_point.shape[1] == 2
+                    point_numbers = list(range(len(cv_point[1:,0])))
+                    P1 = plt.scatter(cv_point[1:,0], cv_point[1:,1], s = 2, marker = 'o', c = point_numbers)
+                    P_init = plt.scatter(cv_point[0,0], cv_point[0,1], s = 20, marker = 'x', c = "k")
+            cbar = plt.colorbar(CS)
+            cbar.ax.tick_params(labelsize=8) 
+            cbar.ax.set_title('kcal/mol',fontsize=8)
+            if op_in["plm_out"] is not None:
+                plt.legend([P_init],["init"])
+                cbar2 = plt.colorbar(P1)
+                cbar2.ax.tick_params(labelsize=4) 
+                cbar2.ax.set_title('steps',fontsize=6)
+            plt.xlabel(r'CV index 1')
+            plt.ylabel(r'CV index 2')
+            plt.savefig(mcmc_2cv_fig,dpi=600,bbox_inches='tight')
+            for iii in range(len(fourdata_2cv)):
+                fig = plt.figure(figsize=(8, 6))
+                cmap = plt.cm.get_cmap("jet_r")
+                # Define percentiles for the levels
+                upper_perc = np.percentile(np.unique(fourdata_2cv[iii]), 99)
+                CS = plt.contourf(xedges,yedges,fourdata_2cv[iii].reshape(bins,bins),levels = np.linspace(0,upper_perc,101),cmap=cmap,extend="max")
+                if op_in["plm_out"] is not None:
+                    for cv_output in op_in["plm_out"]:
+                        cv_point = np.loadtxt(cv_output)
+                        assert cv_point.shape[1] == 2
+                        point_numbers = list(range(len(cv_point[1:,0])))
+                        P1 = plt.scatter(cv_point[1:,0], cv_point[1:,1], s = 2, marker = 'o', c = point_numbers)
+                        P_init = plt.scatter(cv_point[0,0], cv_point[0,1], s = 20, marker = 'x', c = "k")
+                cbar = plt.colorbar(CS)
+                cbar.ax.tick_params(labelsize=8) 
+                cbar.ax.set_title('kcal/mol',fontsize=8)
+                if op_in["plm_out"] is not None:
+                    plt.legend([P_init],["init"])
+                    cbar2 = plt.colorbar(P1)
+                    cbar2.ax.tick_params(labelsize=4) 
+                    cbar2.ax.set_title('steps',fontsize=6)
+                plt.xlabel(r'CV index 1')
+                plt.ylabel(r'CV index 2')
+                plt.savefig(mcmc_2cv_fig_separate.format(tag=iii),dpi=600,bbox_inches='tight')
+                
+        op_out = OPIO(
+            {
+               "mcmc_fig": task_path
+            }
+        )
+        return op_out
\ No newline at end of file
diff --git a/rid/op/mcmc_run.py b/rid/op/mcmc_run.py
new file mode 100644
index 00000000..5cf37670
--- /dev/null
+++ b/rid/op/mcmc_run.py
@@ -0,0 +1,195 @@
+from typing import List, Optional, Dict
+from pathlib import Path
+import numpy as np
+from dflow.python import (
+    OP,
+    OPIO,
+    OPIOSign,
+    Artifact,
+    Parameter,
+    BigParameter
+)
+from rid.utils import save_txt, set_directory
+from rid.mcmc.walker import Walker, my_hist1d, my_hist2d, my_hist1d_path, my_hist2d_path
+from rid.select.model_devi import test_ef
+from rid.common.tensorflow.graph import load_graph
+try:
+    import tensorflow.compat.v1 as tf
+    tf.disable_v2_behavior()
+except ImportError:
+    import tensorflow as tf
+from rid.constants import mcmc_1cv_name, mcmc_1cv_dir_name, mcmc_2cv_name, kb, f_cvt, kcal2kj
+import os
+
+class MCMCRun(OP):
+    """
+    `MCMC_Run` performs Markov chain Monte carlo sampling to reduce the dimension of Free energy surface (represented by neural network)
+    """
+
+    @classmethod
+    def get_input_sign(cls):
+        return OPIOSign(
+            {
+                "task_names": str,
+                "mcmc_config": BigParameter(Dict),
+                "models": Artifact(Path, archive = None)
+            }
+        )
+
+    @classmethod
+    def get_output_sign(cls):
+        return OPIOSign(
+            {
+                "mcmc_1cv": Artifact(Path, archive = None),
+                "mcmc_2cv": Artifact(Path, archive = None)
+            }
+        )
+
+    @OP.exec_sign_check
+    def execute(
+        self,
+        op_in: OPIO,
+    ) -> OPIO:
+
+        r"""Execute the OP.
+        
+        Parameters
+        ----------
+        op_in : dict
+            Input dict with components:
+
+            - `task_name`: str,
+            - `models`: Artifact(List[Path], optional=True),
+          
+        Returns
+        -------
+            Output dict with components:
+        
+            - `fes_fig`: (`Artifact(Path)`)
+        """
+        graph = load_graph(op_in["models"])
+        mcmc_config = op_in["mcmc_config"]
+        temperature = mcmc_config["temperature"]
+        fd = mcmc_config["cv_dimension"]
+        ns = mcmc_config["numb_steps"]
+        nw = mcmc_config["numb_walkers"]
+        # kinetic enery in eV
+        kbT = kb * temperature
+        beta = 1.0 / kbT
+        if "cv_upper_bound" in mcmc_config:
+            cv_upper = mcmc_config["cv_upper_bound"]
+        else:
+            cv_upper = 2*np.pi
+        if "cv_lower_bound" in mcmc_config:
+            cv_lower = mcmc_config["cv_lower_bound"]
+        else:
+            cv_lower = 0
+        proj_info = mcmc_config["proj_info"]
+        proj_mode = proj_info["proj_mode"]
+        proj_cv_index = proj_info["proj_cv_index"]
+        if "path_list" in proj_info:
+            path_list = np.array(proj_info["path_list"])
+            path_lm = proj_info["path_lm"]
+        cv_type = mcmc_config["cv_type"]
+        bins = mcmc_config["bins"]
+        
+        if cv_type == "dih":
+            xx = np.linspace(0,2* np.pi, bins)
+            yy = np.linspace(0,2* np.pi, bins)
+            pp_hist = np.zeros((fd, len(xx)))
+            pp_hist2d = np.zeros((1, len(xx), len(yy)))
+            delta = 2.0 * np.pi / (bins-1)
+        elif cv_type == "dis":
+            xx = np.linspace(0,10, bins)
+            yy = np.linspace(0,10, bins)
+            if proj_mode == "cv":
+                pp_hist = np.zeros((fd, len(xx)))
+            elif proj_mode == "path":
+                pp_hist = np.zeros((2, len(xx)))
+            pp_hist2d = np.zeros((1, len(xx), len(yy)))
+            delta = 10.0 / (bins-1)
+        else:
+            raise ValueError("Undefined cv type, only support 'dih' and 'dis' type")
+        
+        task_path = Path(op_in["task_names"])
+        task_path.mkdir(exist_ok=True, parents=True)
+
+        mcmc_2cv_path = None
+        with set_directory(task_path):
+            with tf.Session(graph = graph) as sess:        
+                walker = Walker(fd, nw, sess, cv_type, cv_lower=cv_lower, cv_upper=cv_upper)
+                for _ in range(10000):
+                    pp, ee, ff = walker.sample(test_ef)
+
+                for ii in range(ns+1):
+                    pp, ee, ff = walker.sample(test_ef)
+                        
+                    if proj_mode == "cv":
+                        # project on 1D CV
+                        pp_hist_new = my_hist1d(pp, xx, delta, fd)
+                        pp_hist = (pp_hist * ii + pp_hist_new) / (ii+1)
+                        if not os.path.exists(mcmc_1cv_dir_name):
+                            os.makedirs(mcmc_1cv_dir_name)
+                        if np.mod(ii,int(ns/5)) == 0:
+                            zz = -np.log(pp_hist+1e-7)/beta
+                            # convert ev to kcal/mol
+                            zz *= f_cvt/kcal2kj
+                            zz = zz - np.min(zz)      
+                            for jj in range(fd):
+                                fp = open(mcmc_1cv_dir_name+"/"+mcmc_1cv_name.format(tag=jj), "a")
+                                for temp in zz[jj]:
+                                    fp.write(str(temp)+'    ')
+                                fp.write('\n')
+                                fp.close()
+                        # project on 2D CV
+                        assert len(proj_cv_index) == 2
+                        cv1 = proj_cv_index[0]
+                        cv2 = proj_cv_index[1]
+                        ##certain 2d
+                        pp_hist_new2d = my_hist2d(pp, xx, yy, delta, cv1, cv2)
+                        pp_hist2d = (pp_hist2d * ii + pp_hist_new2d) / (ii+1)
+                        if ii == ns:
+                            zz2d = np.transpose(-np.log(pp_hist2d+1e-10), (0,2,1))/beta
+                            # convert ev to kcal/mol
+                            zz2d *= f_cvt/kcal2kj
+                            zz2d = zz2d - np.min(zz2d)
+                            np.savetxt(mcmc_2cv_name,zz2d[0])
+                    elif proj_mode == "path":
+                        # project on 1D CV
+                        pp_hist_new = my_hist1d_path(pp, xx, delta, path_lm, path_list, proj_cv_index)
+                        pp_hist = (pp_hist * ii + pp_hist_new) / (ii+1)
+                        if not os.path.exists(mcmc_1cv_dir_name):
+                            os.makedirs(mcmc_1cv_dir_name)
+                        if np.mod(ii,int(ns/5)) == 0:
+                            zz = -np.log(pp_hist+1e-7)/beta
+                            # convert ev to kcal/mol
+                            zz *= f_cvt/kcal2kj
+                            zz = zz - np.min(zz)    
+                            # iterate over 2 path CV  
+                            for jj in range(2):
+                                fp = open(mcmc_1cv_dir_name+"/"+mcmc_1cv_name.format(tag=jj), "a")
+                                for temp in zz[jj]:
+                                    fp.write(str(temp)+'    ')
+                                fp.write('\n')
+                                fp.close()
+                        # project on 2D CV
+                        pp_hist_new2d_path = my_hist2d_path(pp, xx, yy, delta, path_lm, path_list, proj_cv_index)
+                        pp_hist2d = (pp_hist2d * ii + pp_hist_new2d_path) / (ii+1)
+                        if ii == ns:
+                            zz2d = np.transpose(-np.log(pp_hist2d+1e-10), (0,2,1))/beta
+                            # convert ev to kcal/mol
+                            zz2d *= f_cvt/kcal2kj
+                            zz2d = zz2d - np.min(zz2d)
+                            np.savetxt(mcmc_2cv_name,zz2d[0])
+                            
+        if os.path.exists(task_path.joinpath(mcmc_2cv_name)):
+            mcmc_2cv_path = task_path.joinpath(mcmc_2cv_name)
+            
+            
+        op_out = OPIO(
+            {
+               "mcmc_1cv": task_path.joinpath(mcmc_1cv_dir_name),
+               "mcmc_2cv": mcmc_2cv_path
+            }
+        )
+        return op_out
\ No newline at end of file
diff --git a/rid/op/prep_data.py b/rid/op/prep_data.py
index 80c7d7f0..9b314764 100644
--- a/rid/op/prep_data.py
+++ b/rid/op/prep_data.py
@@ -43,10 +43,11 @@ def execute(
         ) -> OPIO:
         cv_forces = []
         for idx in range(len(op_in["cv_forces"])):
-            cv_force = load_txt(op_in["cv_forces"][idx])
-            cv_forces = np.append(cv_forces, cv_force)
+            if op_in["cv_forces"][idx]:
+                cv_force = load_txt(op_in["cv_forces"][idx])
+                cv_forces = np.append(cv_forces, cv_force)
         if op_in["cv_forces"]:
-            cv_forces = np.reshape(cv_forces, [len(op_in["cv_forces"]), -1])
+            cv_forces = np.reshape(cv_forces, [-1, len(cv_force)])
             data = cv_forces
         else:
             data = np.array([])
diff --git a/rid/op/prep_exploration.py b/rid/op/prep_exploration.py
index a749614c..f2f9c0f9 100644
--- a/rid/op/prep_exploration.py
+++ b/rid/op/prep_exploration.py
@@ -3,7 +3,8 @@
     OP,
     OPIO,
     OPIOSign,
-    Artifact
+    Artifact,
+    BigParameter
 )
 
 from typing import List, Dict, Union
@@ -51,8 +52,8 @@ def get_input_sign(cls):
                 "cv_file": Artifact(List[Path], optional=True),
                 "trust_lvl_1": float,
                 "trust_lvl_2": float,
-                "exploration_config": Dict,
-                "cv_config": Dict,
+                "exploration_config": BigParameter(Dict),
+                "cv_config": BigParameter(Dict),
                 "task_name": str,
                 "block_tag": str
             }
@@ -62,7 +63,7 @@ def get_input_sign(cls):
     def get_output_sign(cls):
         return OPIOSign(
             {
-                "task_path": Artifact(Path),
+                "task_path": Artifact(Path, archive = None),
                 "cv_dim": int
             }
         )
diff --git a/rid/op/prep_label.py b/rid/op/prep_label.py
index 33f52f41..dfc99b9c 100644
--- a/rid/op/prep_label.py
+++ b/rid/op/prep_label.py
@@ -3,7 +3,8 @@
     OPIO,
     OPIOSign,
     Artifact,
-    Parameter
+    Parameter,
+    BigParameter
 )
 
 import json
@@ -29,7 +30,7 @@ def get_input_sign(cls):
         return OPIOSign(
             {
                 "confs": Artifact(List[Path], optional=True),
-                "conf_tags": Parameter(List, default=[])
+                "conf_tags": Artifact(List[Path], optional=True)
             }
         )
 
@@ -38,7 +39,7 @@ def get_output_sign(cls):
         return OPIOSign(
             {
                 "if_continue": int,
-                "conf_tags": List
+                "conf_tags": BigParameter(List)
             }
         )
 
@@ -69,10 +70,9 @@ def execute(
 
             tags = {}
             for tag in op_in["conf_tags"]:
-                if isinstance(tag, Dict):
-                    tags.update(tag)
-                elif isinstance(tag, str):
-                    tags.update(json.loads(tag))
+                if isinstance(tag,Path):
+                    with open(tag,"r") as f:
+                        tags.update(json.load(f))
                 else:
                     raise RuntimeError("Unkown Error.")
             
@@ -104,9 +104,9 @@ def get_input_sign(cls):
                 "topology": Artifact(Path, optional=True),
                 "conf": Artifact(Path),
                 "cv_file": Artifact(List[Path], optional=True),
-                "label_config": Dict,
-                "cv_config": Dict,
-                "task_name": str,
+                "label_config": BigParameter(Dict),
+                "cv_config": BigParameter(Dict),
+                "task_name": BigParameter(str),
                 "at": Artifact(Path, optional=True)
             }
         )
@@ -115,7 +115,7 @@ def get_input_sign(cls):
     def get_output_sign(cls):
         return OPIOSign(
             {
-                "task_path": Artifact(Path),
+                "task_path": Artifact(Path, archive = None),
             }
         )
 
diff --git a/rid/op/prep_rid.py b/rid/op/prep_rid.py
index 320d7b2a..1eae2d18 100644
--- a/rid/op/prep_rid.py
+++ b/rid/op/prep_rid.py
@@ -6,10 +6,11 @@
     OP,
     OPIO,
     OPIOSign,
-    Artifact
+    Artifact,
+    BigParameter
 )
 from rid.utils import load_json
-from rid.constants import model_tag_fmt, init_conf_name, init_input_name, walker_tag_fmt
+from rid.constants import model_tag_fmt, init_conf_gmx_name, init_conf_lmp_name,init_input_name, walker_tag_fmt
 
 
 logging.basicConfig(
@@ -21,22 +22,39 @@
 logger = logging.getLogger(__name__)
 
 
-def prep_confs(confs, numb_walkers):
+def prep_confs(confs, numb_walkers,sampler_type):
     numb_confs = len(confs)
     conf_list = []
-    if numb_confs < numb_walkers:
-        logger.info("Number of confs is smaller than number of walkers. Copy replicas up to number of walkers.")
+    if sampler_type == "gmx":
+        if numb_confs < numb_walkers:
+            logger.info("Number of confs is smaller than number of walkers. Copy replicas up to number of walkers.")
+            for idx in range(numb_walkers):
+                shutil.copyfile(confs[idx%numb_confs], init_conf_gmx_name.format(idx=idx))
+        elif numb_confs > numb_walkers:
+            logger.info("Number of confs is greater than number of walkers. Only use the fist `numb_walkers` confs.")
+            for idx in range(numb_walkers):
+                shutil.copyfile(confs[idx], init_conf_gmx_name.format(idx=idx))
+        else:
+            for idx in range(numb_walkers):
+                shutil.copyfile(confs[idx], init_conf_gmx_name.format(idx=idx))
         for idx in range(numb_walkers):
-            shutil.copyfile(confs[idx%numb_confs], init_conf_name.format(idx=idx))
-    elif numb_confs > numb_walkers:
-        logger.info("Number of confs is greater than number of walkers. Only use the fist `numb_walkers` confs.")
+            conf_list.append(Path(init_conf_gmx_name.format(idx=idx)))
+    elif sampler_type == "lmp":
+        if numb_confs < numb_walkers:
+            logger.info("Number of confs is smaller than number of walkers. Copy replicas up to number of walkers.")
+            for idx in range(numb_walkers):
+                shutil.copyfile(confs[idx%numb_confs], init_conf_lmp_name.format(idx=idx))
+        elif numb_confs > numb_walkers:
+            logger.info("Number of confs is greater than number of walkers. Only use the fist `numb_walkers` confs.")
+            for idx in range(numb_walkers):
+                shutil.copyfile(confs[idx], init_conf_lmp_name.format(idx=idx))
+        else:
+            for idx in range(numb_walkers):
+                shutil.copyfile(confs[idx], init_conf_lmp_name.format(idx=idx))
         for idx in range(numb_walkers):
-            shutil.copyfile(confs[idx], init_conf_name.format(idx=idx))
+            conf_list.append(Path(init_conf_lmp_name.format(idx=idx)))
     else:
-        for idx in range(numb_walkers):
-            shutil.copyfile(confs[idx], init_conf_name.format(idx=idx))
-    for idx in range(numb_walkers):
-        conf_list.append(Path(init_conf_name.format(idx=idx)))
+        raise ValueError("Invalid sampler type, only support gmx and lmp")
     return conf_list
 
 
@@ -65,12 +83,12 @@ def get_output_sign(cls):
                 "numb_iters": int,
                 "numb_walkers": int,
                 "numb_models": int,
-                "confs": Artifact(List[Path]),
+                "confs": Artifact(List[Path],archive = None),
                 "walker_tags": List,
                 "model_tags": List,
                 
-                "exploration_config": Dict,
-                "cv_config": Dict,
+                "exploration_config": BigParameter(Dict),
+                "cv_config": BigParameter(Dict),
                 "trust_lvl_1": List[float],
                 "trust_lvl_2": List[float],
                 "cluster_threshold": List[float],
@@ -80,11 +98,13 @@ def get_output_sign(cls):
                 "numb_cluster_lower": int,
                 "max_selection": int,
                 "numb_cluster_threshold": int,
+                "std_threshold": float,
                 "dt": float,
                 "output_freq": float,
                 "slice_mode": str,
-                "label_config": Dict,
-                "train_config": Dict
+                "type_map": List,
+                "label_config": BigParameter(Dict),
+                "train_config": BigParameter(Dict)
             }
         )
 
@@ -142,7 +162,10 @@ def execute(
         train_config = jdata.pop("Train")
         numb_models = train_config.pop("numb_models")
         numb_iters = jdata.pop("numb_iters")
-        conf_list = prep_confs(op_in["confs"], numb_walkers)
+        exploration_config = jdata.pop("ExploreMDConfig")
+        
+        sampler_type = exploration_config["type"]
+        conf_list = prep_confs(op_in["confs"], numb_walkers, sampler_type)
 
         walker_tags = []
         model_tags = []
@@ -151,7 +174,6 @@ def execute(
         for idx in range(numb_models):
             model_tags.append(model_tag_fmt.format(idx=idx))
         
-        exploration_config = jdata.pop("ExploreMDConfig")
         dt = exploration_config["dt"]
         output_freq = exploration_config["output_freq"]
         cv_config = jdata.pop("CV")
@@ -170,6 +192,13 @@ def execute(
         cluster_threshold = selection_config.pop("cluster_threshold")
         cluster_threshold_list = [cluster_threshold for _ in range(numb_walkers)]
         
+        std_threshold = label_config["std_threshold"]
+        
+        if "type_map" in selection_config:
+            type_map = selection_config["type_map"]
+        else:
+            type_map = []
+        
         op_out = OPIO(
             {
                 "numb_iters": numb_iters,
@@ -190,9 +219,11 @@ def execute(
                 "numb_cluster_lower": selection_config.pop("numb_cluster_lower"),
                 "max_selection": selection_config.pop("max_selection"),
                 "numb_cluster_threshold": selection_config.pop("numb_cluster_threshold"),
+                "std_threshold": std_threshold,
                 "dt": dt,
                 "output_freq": output_freq,
                 "slice_mode": selection_config.pop("slice_mode"),
+                "type_map": type_map,
                 "label_config": label_config,
                 "train_config": train_config
             }
diff --git a/rid/op/prep_select.py b/rid/op/prep_select.py
index 88c058bd..8a9f4ed1 100644
--- a/rid/op/prep_select.py
+++ b/rid/op/prep_select.py
@@ -53,9 +53,9 @@ def get_output_sign(cls):
             {
                 "numb_cluster": int,
                 "cluster_threshold": float,
-                "cluster_fig": Artifact(Path),
-                "cluster_selection_index": Artifact(Path),
-                "cluster_selection_data": Artifact(Path)
+                "cluster_fig": Artifact(Path, archive = None),
+                "cluster_selection_index": Artifact(Path, archive = None),
+                "cluster_selection_data": Artifact(Path, archive = None)
             }
         )
 
diff --git a/rid/op/run_exploration.py b/rid/op/run_exploration.py
index 65df3a07..0193b660 100644
--- a/rid/op/run_exploration.py
+++ b/rid/op/run_exploration.py
@@ -8,7 +8,8 @@
     OP,
     OPIO,
     OPIOSign,
-    Artifact
+    Artifact,
+    BigParameter
 )
 from rid.constants import (
         gmx_conf_name,
@@ -24,12 +25,17 @@
         lmp_conf_out,
         bias_fig,
         model_devi_fig,
+        dp_model_devi_fig,
+        dp_sel_ndx,
+        sel_lmp_name,
         lmp_input_name
     )
 from rid.utils import run_command, set_directory, list_to_string
 from rid.common.lammps.command import final_dump
 from rid.common.sampler.command import get_grompp_cmd, get_mdrun_cmd
+from rid.common.mol_dpdata import slice_dump
 from rid.select.model_devi import make_std
+from rid.mcmc.walker import s,z
 
 
 logging.basicConfig(
@@ -55,7 +61,7 @@ def get_input_sign(cls):
             {
                 "task_path": Artifact(Path),
                 "forcefield": Artifact(Path, optional=True),
-                "exploration_config": Dict,
+                "exploration_config": BigParameter(Dict),
                 "models": Artifact(List[Path], optional=True),
                 "index_file": Artifact(Path, optional=True),
                 "dp_files": Artifact(List[Path], optional=True),
@@ -68,12 +74,17 @@ def get_input_sign(cls):
     def get_output_sign(cls):
         return OPIOSign(
             {
-                "plm_out": Artifact(Path),
-                "bias_fig": Artifact(Path, optional=True),
-                "model_devi_fig": Artifact(Path,optional=True),
-                "md_log": Artifact(Path),
-                "trajectory": Artifact(Path),
-                "conf_out": Artifact(Path)
+                "plm_out": Artifact(Path, archive = None),
+                "bias_fig": Artifact(Path, optional=True, archive = None),
+                "model_devi_fig": Artifact(Path,optional=True, archive = None),
+                "dp_model_devi_fig": Artifact(Path,optional=True, archive = None),
+                "dp_model_devi": Artifact(Path,optional=True, archive = None),
+                "dp_selected_indices": Artifact(Path,optional=True, archive = None),
+                "dp_selected_confs": Artifact(List[Path], archive = None),
+                "projected_fig": Artifact(Path,optional=True, archive = None),
+                "md_log": Artifact(Path, archive = None),
+                "trajectory": Artifact(Path, archive = None),
+                "conf_out": Artifact(Path, archive = None)
             }
         )
 
@@ -159,24 +170,28 @@ def execute(
                 for file in op_in["dp_files"]:
                     if not os.path.islink(file.name):
                         os.symlink(file, file.name)
-            if op_in["index_file"] is not None:
-                if not os.path.islink(op_in["index_file"].name):
-                    os.symlink(op_in["index_file"], op_in["index_file"].name)
+
             if op_in["cv_file"] is not None:
                 for file in op_in["cv_file"]:
                     if file.name != "colvar":
                         if not os.path.islink(file.name):
                             os.symlink(file, file.name)
-            if op_in["inputfile"] is not None:
-                for file in op_in["inputfile"]:
-                    if file.name == inputfile_name:
-                        if not os.path.islink(file.name):
-                            os.symlink(file, file.name)
+                            
+            if op_in["exploration_config"]["type"] == "lmp":
+                if op_in["inputfile"] is not None:
+                    for file in op_in["inputfile"]:
+                        if file.name == inputfile_name:
+                            if not os.path.islink(file.name):
+                                os.symlink(file, file.name)
             
-            if op_in["dp_files"] is not None:
-                for dp_file in op_in["dp_files"]:
-                    if (dp_file.name).endswith("json"):
-                        os.environ["GMX_DEEPMD_INPUT_JSON"] = dp_file.name
+            if op_in["exploration_config"]["type"] == "gmx":
+                if op_in["index_file"] is not None:
+                    if not os.path.islink(op_in["index_file"].name):
+                        os.symlink(op_in["index_file"], op_in["index_file"].name)
+                if op_in["dp_files"] is not None:
+                    for dp_file in op_in["dp_files"]:
+                        if (dp_file.name).endswith("json"):
+                            os.environ["GMX_DEEPMD_INPUT_JSON"] = dp_file.name
                         
             if grompp_cmd is not None:
                 logger.info(list_to_string(grompp_cmd, " "))
@@ -193,6 +208,77 @@ def execute(
             
             bias_fig_file = None
             model_devi_fig_file = None
+            
+            dp_model_devi_fig_file = None
+            dp_model_devi_file = None
+            dp_sel_ndx_file = None
+            dp_conf_list = []
+            if "dp_model_devi_out" in op_in["exploration_config"]:
+                # plot DP potential model deviation
+                dp_model_devi_file_name = op_in["exploration_config"]["dp_model_devi_out"]
+                data = np.loadtxt(dp_model_devi_file_name)
+                f_list = data[:,4]
+                cv_data = np.loadtxt("plm.out")[:,2:]
+                assert (len(cv_data) == len(f_list))
+                trust_lvl_hi = op_in["exploration_config"]["dp_model_devi_f_hi"]
+                trust_lvl_lo = op_in["exploration_config"]["dp_model_devi_f_lo"]
+                max_selection = op_in["exploration_config"]["dp_model_max_selection"]
+                high_devi = np.where((f_list>trust_lvl_lo) & (f_list < trust_lvl_hi))
+                
+                if len(high_devi[0]) > max_selection:
+                    high_devi = np.random.choice(high_devi[0],max_selection,replace=False)
+                else:
+                    high_devi = high_devi[0]
+                np.savetxt(dp_sel_ndx,high_devi,fmt="%d")
+                dp_sel_ndx_file = op_in["task_path"].joinpath(dp_sel_ndx)
+                walker_idx = int(op_in["task_path"].name)
+                slice_dump(dump="out.dump",walker_idx = walker_idx,selected_idx=high_devi, output=sel_lmp_name,style="dpdata")
+                for ii, sel in enumerate(high_devi):
+                    dp_conf_list.append(op_in["task_path"].joinpath(sel_lmp_name.format(walker=walker_idx,idx=sel)))
+                
+                x_list = np.array([i for i in range(len(f_list))])
+                nframes = len(x_list)
+                groups = ["other"]*nframes
+                y_list = [int(j) for j in high_devi]
+                for j in y_list:
+                    groups[j] = "choice"
+                cdict = {"other": 'black', "choice": 'red'}
+                groups = np.array(groups)
+                fig, ax = plt.subplots()
+                for g in ["other","choice"]:
+                    ix = np.where(groups == g)
+                    ax.scatter(x_list[ix], f_list[ix], c = cdict[g], label = g, s = 20)
+                plt.title("max model devi of DP model")
+                plt.legend()
+                plt.xlabel("simu frames")
+                plt.ylabel("max model devi")
+                plt.savefig(dp_model_devi_fig)
+                dp_model_devi_file = op_in["task_path"].joinpath(dp_model_devi_file_name)
+                dp_model_devi_fig_file = op_in["task_path"].joinpath(dp_model_devi_fig)
+            
+            projected_fig_file = None
+            if "proj_info" in op_in["exploration_config"]:
+                proj_info = op_in["exploration_config"]["proj_info"]
+                # plot projected cv fig
+                proj_mode = proj_info["proj_mode"]
+                if proj_mode == "path":
+                    proj_cv_index = proj_info["proj_cv_index"]
+                    proj_lm = proj_info["path_lm"]
+                    path_list = proj_info["path_list"]
+                    cv_data = np.loadtxt("plm.out")[:,2:]
+                    proj_cv_data = cv_data[:,proj_cv_index]                                      
+                    slist = s(proj_cv_data, proj_lm, path_list)
+                    zlist = z(proj_cv_data, proj_lm, path_list)
+                    point_numbers = list(range(len(slist)))
+                    plt.figure(figsize=(10, 8), dpi=100)
+                    plt.scatter(slist,zlist, c = point_numbers)
+                    plt.xlabel("s")
+                    plt.ylabel("z")
+                    plt.title("Trajectory projection")
+                    plt.colorbar(label="steps")
+                    plt.savefig("path.png")
+                    projected_fig_file = op_in["task_path"].joinpath("path.png")
+            
             if op_in["models"] is not None:
                 # plot bias during simulation
                 bias = np.loadtxt(plumed_output_name)[:,1]
@@ -234,6 +320,11 @@ def execute(
                 "plm_out": op_in["task_path"].joinpath(plumed_output_name),
                 "bias_fig": bias_fig_file,
                 "model_devi_fig": model_devi_fig_file,
+                "dp_model_devi_fig": dp_model_devi_fig_file,
+                "dp_model_devi": dp_model_devi_file,
+                "dp_selected_indices": dp_sel_ndx_file,
+                "dp_selected_confs": dp_conf_list,
+                "projected_fig": projected_fig_file,
                 "md_log": op_in["task_path"].joinpath(mdrun_log),
                 "trajectory": op_in["task_path"].joinpath(traj_name),
                 "conf_out": op_in["task_path"].joinpath(conf_out),
diff --git a/rid/op/run_label.py b/rid/op/run_label.py
index a1edc951..ca5627ae 100644
--- a/rid/op/run_label.py
+++ b/rid/op/run_label.py
@@ -7,7 +7,8 @@
     OPIO,
     OPIOSign,
     Artifact,
-    Parameter
+    Parameter,
+    BigParameter
 )
 from rid.constants import (
         gmx_conf_name,
@@ -54,10 +55,10 @@ def get_input_sign(cls):
             {
                 "task_path": Artifact(Path),
                 "forcefield": Artifact(Path, optional=True),
-                "label_config": Dict,
-                "cv_config": Dict,
+                "label_config": BigParameter(Dict),
+                "cv_config": BigParameter(Dict),
                 "at": Artifact(Path, optional=True),
-                "task_name": str,
+                "task_name": BigParameter(str),
                 "tail": Parameter(float, default=0.9),
                 "index_file": Artifact(Path, optional=True),
                 "dp_files": Artifact(List[Path], optional=True),
@@ -70,11 +71,12 @@ def get_input_sign(cls):
     def get_output_sign(cls):
         return OPIOSign(
             {
-                "plm_out": Artifact(Path),
-                "cv_forces": Artifact(Path),
-                "mf_info": Artifact(Path,optional=True),
-                "mf_fig": Artifact(Path),
-                "md_log": Artifact(Path)
+                "plm_out": Artifact(Path, archive = None),
+                "trajectory": Artifact(Path, archive = None),
+                "cv_forces": Artifact(Path, archive = None),
+                "mf_info": Artifact(Path,optional=True, archive = None),
+                "mf_fig": Artifact(Path, archive = None),
+                "md_log": Artifact(Path, archive = None)
             }
         )
 
@@ -181,15 +183,17 @@ def execute(
         traj_path = None
         if op_in["label_config"]["type"] == "gmx":
             mdrun_log = gmx_mdrun_log
-            if os.path.exists(op_in["task_path"].joinpath(gmx_xtc_name)):
-                traj_path = op_in["task_path"].joinpath(gmx_xtc_name)
+            if "traj_output" in op_in["label_config"] and op_in["label_config"]["traj_output"] == "True":
+                if os.path.exists(op_in["task_path"].joinpath(gmx_xtc_name)):
+                    traj_path = op_in["task_path"].joinpath(gmx_xtc_name)
             if op_in["label_config"]["method"] == "constrained":
                 frame_coords = op_in["task_path"].joinpath(gmx_coord_name)
                 frame_forces = op_in["task_path"].joinpath(gmx_force_name)
         elif op_in["label_config"]["type"] == "lmp":
             mdrun_log = lmp_mdrun_log
-            if os.path.exists(op_in["task_path"].joinpath("out.dump")):
-                traj_path = op_in["task_path"].joinpath("out.dump")
+            if "traj_output" in op_in["label_config"] and op_in["label_config"]["traj_output"] == "True":
+                if os.path.exists(op_in["task_path"].joinpath("out.dump")):
+                    traj_path = op_in["task_path"].joinpath("out.dump")
         
         conf = None
         topology = None
@@ -218,6 +222,7 @@ def execute(
         op_out = OPIO(
             {
                 "plm_out": plm_out,
+                "trajectory": traj_path,
                 "cv_forces": mf_out["cv_forces"],
                 "mf_info":  mf_out["mf_info"],
                 "mf_fig": Path(op_in["task_name"]).joinpath(mf_fig),
diff --git a/rid/op/run_model_devi.py b/rid/op/run_model_devi.py
new file mode 100644
index 00000000..3531d0c3
--- /dev/null
+++ b/rid/op/run_model_devi.py
@@ -0,0 +1,123 @@
+from typing import List, Optional, Dict
+from pathlib import Path
+import numpy as np
+from dflow.python import (
+    OP,
+    OPIO,
+    OPIOSign,
+    Artifact,
+    Parameter,
+    BigParameter
+)
+from rid.utils import save_txt, set_directory
+from rid.constants import sel_gro_name, cv_init_label, model_devi_name, model_devi_precision, sel_ndx_name, new_model_devi_fig
+from rid.select.conf_select import select_from_devi
+from rid.common.mol import get_distance_from_atomid
+from rid.select.model_devi import make_std
+import os
+from matplotlib import pyplot as plt
+
+class RunModelDevi(OP):
+
+    """
+    `RunModelDevi` calculates model deviations for each chosen representive cluster frames from `PrepSelect` and select 
+    ones with high uncertainty from them.
+    """
+
+    @classmethod
+    def get_input_sign(cls):
+        return OPIOSign(
+            {
+                "models": Artifact(List[Path]),
+                "plm_out": Artifact(Path),
+                "selected_indices": Artifact(Path, optional=True),
+                "trust_lvl_1": float,
+                "task_name": str,
+                "block_tag": str
+            }
+        )
+
+    @classmethod
+    def get_output_sign(cls):
+        return OPIOSign(
+            {
+                "model_devi": Artifact(Path, optional=True, archive = None),
+                "model_devi_fig": Artifact(Path, optional=True, archive = None)
+            }
+        )
+
+    @OP.exec_sign_check
+    def execute(
+        self,
+        op_in: OPIO,
+    ) -> OPIO:
+
+        r"""Execute the OP.
+        
+        Parameters
+        ----------
+        op_in : dict
+            Input dict with components:
+
+            - `models`: Artifact(List[Path], optional=True),
+            - `plm_out`: Artifact(Path)
+          
+        Returns
+        -------
+            Output dict with components:
+        
+            - `model_devi`: (`Artifact(Path)`) Model deviation values of selected conformation files (`selected_confs`) .
+        """
+        task_path = Path(op_in["task_name"])
+        task_path.mkdir(exist_ok=True, parents=True)
+
+        with set_directory(task_path):
+            # if "units" in cv_config:
+            #     length_units = cv_config["units"]
+            # else:
+            #     length_units = None
+            # if length_units == None or length_units == "nm":
+            #     cv_data = np.loadtxt(op_in["plm_out"])[:,2:]
+            # elif length_units == "A" or length_units == "Angstrom":
+            #     cv_data = 10*np.loadtxt(op_in["plm_out"])[:,2:]
+            # else:
+            #     raise ValueError("Valid length units must be nm or A")
+            trust_lvl = op_in["trust_lvl_1"]
+            cv_data = np.loadtxt(op_in["plm_out"])[:,2:]
+            stds = make_std(cv_data, models=op_in["models"])
+            save_txt("cls_"+model_devi_name, stds, fmt=model_devi_precision)
+            
+            _selected_idx = np.load(op_in["selected_indices"])
+     
+            nframes = cv_data.shape[0]
+            groups = ["other"]*nframes
+            xlist = np.linspace(0,nframes-1,nframes,dtype=int)
+            ylist = [int(j) for j in _selected_idx]                
+            for y_ in ylist:
+                groups[y_] = "train"
+            cdict = {"other": 'black',"train": 'red'}
+            fig, ax = plt.subplots()
+            groups = np.array(groups)
+            
+            for g in ["other","train"]:
+                ix = np.where(groups == g)
+                ax.scatter(xlist[ix], stds[ix], c = cdict[g], label = g, s = 20)
+            plt.axhline(y=trust_lvl, color='b', linestyle='-',label="trust level 1")
+            plt.xlabel("configuration list")
+            plt.ylabel("force std (KJ/(mol*length))")
+            plt.legend()
+            plt.title("prediction error")
+            fig.set_size_inches(8, 6, forward=True)
+            plt.savefig(new_model_devi_fig)
+                
+        model_devi_fig_path = None
+        if os.path.exists(task_path.joinpath(new_model_devi_fig)):
+            model_devi_fig_path = task_path.joinpath(new_model_devi_fig)
+            
+        op_out = OPIO(
+            {
+               "model_devi": task_path.joinpath("cls_"+model_devi_name),
+               "model_devi_fig": model_devi_fig_path
+            }
+        )
+        return op_out
\ No newline at end of file
diff --git a/rid/op/run_select.py b/rid/op/run_select.py
index 79eacf8b..24ab1e91 100644
--- a/rid/op/run_select.py
+++ b/rid/op/run_select.py
@@ -6,13 +6,16 @@
     OPIO,
     OPIOSign,
     Artifact,
-    Parameter
+    Parameter,
+    BigParameter
 )
 from rid.utils import save_txt, set_directory
 from rid.constants import sel_gro_name, sel_lmp_name, cv_init_label, model_devi_name, model_devi_precision, sel_ndx_name
 from rid.select.conf_select import select_from_devi
-from rid.common.mol import slice_xtc, slice_dump
+from rid.common.mol import slice_xtc
+from rid.common.mol_dpdata import slice_dump
 from rid.select.model_devi import make_std
+import json
 
 
 class RunSelect(OP):
@@ -39,6 +42,8 @@ def get_input_sign(cls):
                 "trust_lvl_2": float,
                 "xtc_traj": Artifact(Path),
                 "topology": Artifact(Path),
+                "label_config": BigParameter(Dict),
+                "type_map": Parameter(Optional[List], default=[]),
                 "dt": Parameter(Optional[float], default=None),
                 "output_freq": Parameter(Optional[float], default=None),
                 "slice_mode": Parameter(str, default="gmx")
@@ -49,11 +54,11 @@ def get_input_sign(cls):
     def get_output_sign(cls):
         return OPIOSign(
             {
-                "selected_confs": Artifact(List[Path]),
-                "selected_cv_init": Artifact(List[Path]),
-                "model_devi": Artifact(Path, optional=True),
-                "selected_indices": Artifact(Path),
-                "selected_conf_tags": Dict
+                "selected_confs": Artifact(List[Path], archive = None),
+                "selected_cv_init": Artifact(List[Path], archive = None),
+                "model_devi": Artifact(Path, optional=True, archive = None),
+                "selected_indices": Artifact(Path, archive = None),
+                "selected_conf_tags": Artifact(Path, archive= None)
             }
         )
 
@@ -107,7 +112,7 @@ def execute(
                 save_txt("cls_"+model_devi_name, stds, fmt=model_devi_precision)
                 _selected_idx = select_from_devi(stds, op_in["trust_lvl_1"])
             sel_idx = cls_sel_idx[_selected_idx]
-            save_txt(sel_ndx_name, sel_idx, fmt="%d")
+            np.save(sel_ndx_name, sel_idx)
             sel_data = cls_sel_data[_selected_idx]
             if op_in["slice_mode"] == "gmx":
                 assert op_in["dt"] is not None, "Please provide time step to slice trajectory."
@@ -119,7 +124,12 @@ def execute(
                 slice_xtc(xtc=op_in["xtc_traj"], top=op_in["topology"],
                         walker_idx = walker_idx, selected_idx=sel_idx, output=sel_gro_name, style="mdtraj")
             elif op_in["slice_mode"] == "dpdata":
-                slice_dump(dump=op_in["xtc_traj"],walker_idx = walker_idx,selected_idx=sel_idx, output=sel_lmp_name, style="dpdata")
+                if op_in["label_config"]["type"] == "gmx":
+                    slice_dump(dump=op_in["xtc_traj"],walker_idx = walker_idx,selected_idx=sel_idx, output=sel_gro_name, type_map = op_in["type_map"], style="dpdata")
+                elif op_in["label_config"]["type"] == "lmp":
+                    slice_dump(dump=op_in["xtc_traj"],walker_idx = walker_idx,selected_idx=sel_idx, output=sel_lmp_name, type_map = op_in["type_map"],style="dpdata")
+                else:
+                    raise ValueError("Invalid labeling type, only support gmx and lmp")
             else:
                 raise RuntimeError("Unknown Style for Slicing Trajectory.")
             conf_list = []
@@ -127,13 +137,19 @@ def execute(
             conf_tags = {}
             for ii, sel in enumerate(sel_idx):
                 if op_in["slice_mode"] == "dpdata":
-                    conf_list.append(task_path.joinpath(sel_lmp_name.format(walker=walker_idx,idx=sel)))
-                    conf_tags[sel_lmp_name.format(walker = walker_idx,idx=sel)] = f"{op_in['task_name']}_{sel}"
+                    if op_in["label_config"]["type"] == "gmx":
+                        conf_list.append(task_path.joinpath(sel_gro_name.format(walker=walker_idx,idx=sel)))
+                        conf_tags[sel_gro_name.format(walker = walker_idx,idx=sel)] = f"{op_in['task_name']}_{sel}"
+                    elif op_in["label_config"]["type"] == "lmp":
+                        conf_list.append(task_path.joinpath(sel_lmp_name.format(walker=walker_idx,idx=sel)))
+                        conf_tags[sel_lmp_name.format(walker = walker_idx,idx=sel)] = f"{op_in['task_name']}_{sel}"
                 elif op_in["slice_mode"] == "gmx" or op_in["slice_mode"] == "mdtraj" :
                     conf_list.append(task_path.joinpath(sel_gro_name.format(walker=walker_idx,idx=sel)))
                     conf_tags[sel_gro_name.format(walker = walker_idx,idx=sel)] = f"{op_in['task_name']}_{sel}"
                 save_txt(cv_init_label.format(walker=walker_idx,idx=sel), sel_data[ii])
                 cv_init_list.append(task_path.joinpath(cv_init_label.format(walker=walker_idx,idx=sel)))
+            with open("conf.json", "w") as f:
+                json.dump(conf_tags,f)
             
         op_out = OPIO(
             {
@@ -141,7 +157,7 @@ def execute(
                "selected_cv_init": cv_init_list,
                "model_devi": task_path.joinpath("cls_"+model_devi_name),
                "selected_indices": task_path.joinpath(sel_ndx_name),
-               "selected_conf_tags": conf_tags
+               "selected_conf_tags": task_path.joinpath("conf.json")
             }
         )
         return op_out
\ No newline at end of file
diff --git a/rid/op/run_train.py b/rid/op/run_train.py
index 0e519150..39140a26 100644
--- a/rid/op/run_train.py
+++ b/rid/op/run_train.py
@@ -6,7 +6,8 @@
     OPIO,
     OPIOSign,
     Artifact,
-    Parameter
+    Parameter,
+    BigParameter
 )
 from rid.constants import tf_model_name, train_fig, train_log
 from rid.nn.train_net import train
@@ -28,7 +29,7 @@ def get_input_sign(cls):
                 "model_tag": str,
                 "angular_mask": List,
                 "data": Artifact(Path),
-                "train_config": Dict
+                "train_config": BigParameter(Dict)
             }
         )
 
diff --git a/rid/select/cluster.py b/rid/select/cluster.py
index e43924c0..c3119e67 100644
--- a/rid/select/cluster.py
+++ b/rid/select/cluster.py
@@ -105,6 +105,26 @@ def mk_cluster(dist, distance_threshold):
     return cluster.labels_
 
 
+def chooseClusterCenter(dist:np.ndarray, conf_ids:list):
+    id_min = conf_ids[0]
+    min_loss = None
+    # n_tot_frame = dist.shape[0]
+    for id in conf_ids:
+        rmsd = dist[id]
+        loss = 0
+        for i, conf in enumerate(conf_ids):
+            rms = rmsd[conf]
+            loss += rms * rms
+        loss = np.sqrt(loss)
+        if min_loss == None:
+            min_loss = loss
+        else:
+            if loss < min_loss:
+                id_min = id
+                min_loss = loss
+    return [id_min]
+
+
 def sel_from_cluster(cvs, threshold, task_path, angular_mask=None, weights=None, max_selection=1000):
     if len(cvs) <= 1:
         return cvs
@@ -133,7 +153,8 @@ def sel_from_cluster(cvs, threshold, task_path, angular_mask=None, weights=None,
     cls_sel = []
     np.random.seed(seed=None)
     for cluster, _ in cls_map:
-        _ret = np.random.choice(cluster, 1)
+        # _ret = np.random.choice(cluster, 1)
+        _ret = chooseClusterCenter(dist, cluster)
         cls_sel.append(_ret[0])
     if len(cls_sel) > max_selection:
         cls_sel = cls_sel[:max_selection]
diff --git a/rid/superop/blocks.py b/rid/superop/blocks.py
index 51ec9bac..6a2709fe 100755
--- a/rid/superop/blocks.py
+++ b/rid/superop/blocks.py
@@ -40,7 +40,9 @@ def __init__(
         label_op: OP,
         data_op: OP,
         train_op: OP,
+        model_devi_op: OP,
         train_config: Dict,
+        model_devi_config: Dict,
         upload_python_package = None,
         retry_times = None
     ):
@@ -59,10 +61,12 @@ def __init__(
             "numb_cluster_upper": InputParameter(type=float),
             "numb_cluster_lower": InputParameter(type=float),
             "max_selection": InputParameter(type=int),
+            "std_threshold": InputParameter(type=float, value=5.0),
             "dt": InputParameter(type=float, value=0.02),
             "output_freq": InputParameter(type=float, value=2500),
             "slice_mode": InputParameter(type=str, value="gmx"),
             "label_config": InputParameter(type=Dict),
+            "type_map": InputParameter(type=List, value = []),
             "tail": InputParameter(type=float, value=0.9),
             "train_config": InputParameter(type=Dict)
         }        
@@ -107,7 +111,9 @@ def __init__(
             label_op,
             data_op,
             train_op,
+            model_devi_op,
             train_config,
+            model_devi_config,
             upload_python_package = upload_python_package,
             retry_times = retry_times
         )            
@@ -140,7 +146,9 @@ def _first_run_block(
         label_op: OP,
         data_op: OP,
         train_op: OP,
+        model_devi_op: OP,
         train_config : Dict,
+        model_devi_config: Dict,
         upload_python_package : str = None,
         retry_times: int = None
     ):
@@ -173,6 +181,7 @@ def _first_run_block(
         "Selection",
         template=select_op,
         parameters={
+            "label_config": block_steps.inputs.parameters["label_config"],
             "trust_lvl_1" : block_steps.inputs.parameters["trust_lvl_1"],
             "trust_lvl_2": block_steps.inputs.parameters["trust_lvl_2"],
             "cluster_threshold": block_steps.inputs.parameters["cluster_threshold"],
@@ -184,6 +193,7 @@ def _first_run_block(
             "dt": block_steps.inputs.parameters["dt"],
             "output_freq": block_steps.inputs.parameters["output_freq"],
             "slice_mode": block_steps.inputs.parameters["slice_mode"],
+            "type_map": block_steps.inputs.parameters["type_map"],
             "if_make_threshold": True,
             "task_names" : block_steps.inputs.parameters['walker_tags'],
             "block_tag" : block_steps.inputs.parameters['block_tag'],
@@ -205,8 +215,8 @@ def _first_run_block(
             "label_config": block_steps.inputs.parameters['label_config'],
             "cv_config": block_steps.inputs.parameters['cv_config'],
             "tail": block_steps.inputs.parameters['tail'],
-            "conf_tags" : selection.outputs.parameters['selected_conf_tags'],
             "block_tag" : block_steps.inputs.parameters['block_tag'],
+            "std_threshold": block_steps.inputs.parameters["std_threshold"]
         },
         artifacts={
             "topology": block_steps.inputs.artifacts["topology"],
@@ -216,7 +226,8 @@ def _first_run_block(
             "index_file": block_steps.inputs.artifacts['index_file'],
             "inputfile": block_steps.inputs.artifacts['inputfile'],
             "dp_files": block_steps.inputs.artifacts['dp_files'],
-            "cv_file": block_steps.inputs.artifacts['cv_file']
+            "cv_file": block_steps.inputs.artifacts['cv_file'],
+            "conf_tags" : selection.outputs.artifacts['selected_conf_tags']
         },
         key = '{}-label'.format(block_steps.inputs.parameters['block_tag'])
     )
@@ -262,6 +273,38 @@ def _first_run_block(
         **train_config,
     )
     block_steps.add(train)
+    
+    model_devi_config = deepcopy(model_devi_config)
+    model_devi_template_config = model_devi_config.pop('template_config')
+    model_devi_executor = init_executor(model_devi_config.pop('executor'))
+    deviation = Step(
+        "ModelDeviation",
+        template=PythonOPTemplate(
+            model_devi_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices("{{item}}",
+                input_parameter=["trust_lvl_1","task_name"],
+                input_artifact=["plm_out","selected_indices"],
+                output_artifact=["model_devi","model_devi_fig"]),
+            **model_devi_template_config,
+        ),
+        parameters={
+            "trust_lvl_1": block_steps.inputs.parameters["trust_lvl_1"],
+            "block_tag" : block_steps.inputs.parameters['block_tag'],
+            "task_name": block_steps.inputs.parameters['walker_tags'],
+        },
+        artifacts={
+            "models" : train.outputs.artifacts["model"],
+            "plm_out": exploration.outputs.artifacts["plm_out"],
+            "selected_indices": selection.outputs.artifacts["selected_indices"]
+        },
+        executor = model_devi_executor,
+        with_param=argo_range(argo_len(block_steps.inputs.parameters["walker_tags"])),
+        key = '{}-model-devi'.format(block_steps.inputs.parameters['block_tag'])+"-{{item}}",
+        **model_devi_config
+    )
+    block_steps.add(deviation)
 
     block_steps.outputs.artifacts["models"]._from = train.outputs.artifacts["model"]
     block_steps.outputs.artifacts["data"]._from = gen_data.outputs.artifacts["data"]
@@ -288,8 +331,10 @@ def __init__(
         data_op: OP,
         adjust_lvl_op: OP,
         train_op: OP,  
+        model_devi_op: OP,
         adjust_lvl_config: Dict,
         train_config: Dict,
+        model_devi_config: Dict,
         upload_python_package = None,
         retry_times = None
     ):
@@ -309,12 +354,14 @@ def __init__(
             "weights": InputParameter(type=Optional[Union[np.ndarray, List]]),
             "max_selection": InputParameter(type=int),
             "numb_cluster_threshold": InputParameter(type=float, value=30),
+            "std_threshold": InputParameter(type=float, value=5.0),
             "dt": InputParameter(type=float, value=0.02),
             "output_freq": InputParameter(type=float, value=2500),
             "slice_mode": InputParameter(type=str, value="gmx"),
             "label_config": InputParameter(type=Dict),
             "tail": InputParameter(type=float, value=0.9),
             "train_config": InputParameter(type=Dict),
+            "type_map": InputParameter(type=List, value=[]),
             "adjust_amplifier": InputParameter(type=float, value=1.5),
             "max_level_multiple": InputParameter(type=float, value=8.0),
         }        
@@ -363,8 +410,10 @@ def __init__(
             data_op,
             adjust_lvl_op,
             train_op,
+            model_devi_op,
             adjust_lvl_config,
             train_config,
+            model_devi_config,
             upload_python_package = upload_python_package,
             retry_times = retry_times
         )            
@@ -398,8 +447,10 @@ def _iter_block(
         data_op: OP,
         adjust_lvl_op: OP,
         train_op: OP,
+        model_devi_op: OP,
         adjust_lvl_config : Dict,
         train_config : Dict,
+        model_devi_config: Dict,
         upload_python_package : str = None,
         retry_times: int = None
     ):
@@ -433,6 +484,7 @@ def _iter_block(
         "Selection",
         template=select_op,
         parameters={
+            "label_config": block_steps.inputs.parameters["label_config"],
             "trust_lvl_1" : block_steps.inputs.parameters["init_trust_lvl_1"],
             "trust_lvl_2": block_steps.inputs.parameters["init_trust_lvl_2"],
             "cluster_threshold": block_steps.inputs.parameters["cluster_threshold"],
@@ -442,6 +494,7 @@ def _iter_block(
             "dt": block_steps.inputs.parameters["dt"],
             "output_freq": block_steps.inputs.parameters["output_freq"],
             "slice_mode": block_steps.inputs.parameters["slice_mode"],
+            "type_map": block_steps.inputs.parameters["type_map"],
             "if_make_threshold": False,
             "task_names" : block_steps.inputs.parameters['walker_tags'],
             "block_tag" : block_steps.inputs.parameters['block_tag'],
@@ -497,8 +550,8 @@ def _iter_block(
             "label_config": block_steps.inputs.parameters['label_config'],
             "cv_config": block_steps.inputs.parameters['cv_config'],
             "tail": block_steps.inputs.parameters['tail'],
-            "conf_tags" : selection.outputs.parameters['selected_conf_tags'],
             "block_tag" : block_steps.inputs.parameters['block_tag'],
+            "std_threshold": block_steps.inputs.parameters["std_threshold"]
         },
         artifacts={
             "topology": block_steps.inputs.artifacts["topology"],
@@ -508,7 +561,8 @@ def _iter_block(
             "at": selection.outputs.artifacts["selected_cv_init"],
             "index_file": block_steps.inputs.artifacts['index_file'],
             "dp_files": block_steps.inputs.artifacts['dp_files'],
-            "cv_file": block_steps.inputs.artifacts['cv_file']
+            "cv_file": block_steps.inputs.artifacts['cv_file'],
+            "conf_tags" : selection.outputs.artifacts['selected_conf_tags']
         },
         key = '{}-label'.format(block_steps.inputs.parameters['block_tag'])
     )
@@ -554,6 +608,38 @@ def _iter_block(
         **train_config,
     )
     block_steps.add(train)
+    
+    model_devi_config = deepcopy(model_devi_config)
+    model_devi_template_config = model_devi_config.pop('template_config')
+    model_devi_executor = init_executor(model_devi_config.pop('executor'))
+    deviation = Step(
+        "ModelDeviation",
+        template=PythonOPTemplate(
+            model_devi_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices("{{item}}",
+                input_parameter=["trust_lvl_1","task_name"],
+                input_artifact=["plm_out","selected_indices"],
+                output_artifact=["model_devi","model_devi_fig"]),
+            **model_devi_template_config,
+        ),
+        parameters={
+            "trust_lvl_1": block_steps.inputs.parameters["trust_lvl_1"],
+            "block_tag" : block_steps.inputs.parameters['block_tag'],
+            "task_name": block_steps.inputs.parameters['walker_tags'],
+        },
+        artifacts={
+            "models" : train.outputs.artifacts["model"],
+            "plm_out": exploration.outputs.artifacts["plm_out"],
+            "selected_indices": selection.outputs.artifacts["selected_indices"]
+        },
+        executor = model_devi_executor,
+        with_param=argo_range(argo_len(block_steps.inputs.parameters["walker_tags"])),
+        key = '{}-model-devi'.format(block_steps.inputs.parameters['block_tag'])+"-{{item}}",
+        **model_devi_config
+    )
+    block_steps.add(deviation)
 
     block_steps.outputs.artifacts["models"]._from = train.outputs.artifacts["model"]
     block_steps.outputs.artifacts["data"]._from = gen_data.outputs.artifacts["data"]
diff --git a/rid/superop/exploration.py b/rid/superop/exploration.py
index 7a85eb5c..7d45071b 100755
--- a/rid/superop/exploration.py
+++ b/rid/superop/exploration.py
@@ -128,7 +128,11 @@ def _exploration(
     run_template_config = run_config.pop('template_config')
     prep_executor = init_executor(prep_config.pop('executor'))
     run_executor = init_executor(run_config.pop('executor'))
-    prep_exploration = Step(
+    prep_merge = False
+    if prep_executor is not None:
+        prep_merge = prep_executor.merge_sliced_step
+    if prep_merge:
+        prep_exploration = Step(
         'prep-exploration',
         template=PythonOPTemplate(
             prep_exploration_op,
@@ -159,11 +163,48 @@ def _exploration(
         key = step_keys["prep_exploration"]+"-{{item}}",
         with_param=argo_range(argo_len(exploration_steps.inputs.parameters['task_names'])),
         executor = prep_executor,
-        **prep_config,
+        **prep_config
+    )
+    else:
+        prep_exploration = Step(
+            'prep-exploration',
+            template=PythonOPTemplate(
+                prep_exploration_op,
+                python_packages = upload_python_package,
+                retry_on_transient_error = retry_times,
+                slices=Slices(sub_path = True,
+                    input_parameter=["task_name", "trust_lvl_1", "trust_lvl_2"],
+                    input_artifact=["conf"],
+                    output_artifact=["task_path"],
+                    output_parameter=["cv_dim"]
+                ),
+                **prep_template_config,
+            ),
+            parameters={
+                "trust_lvl_1" : exploration_steps.inputs.parameters['trust_lvl_1'],
+                "trust_lvl_2": exploration_steps.inputs.parameters['trust_lvl_2'],
+                "exploration_config" : exploration_steps.inputs.parameters['exploration_config'],
+                "cv_config" : exploration_steps.inputs.parameters['cv_config'],
+                "task_name": exploration_steps.inputs.parameters['task_names'],
+                "block_tag": exploration_steps.inputs.parameters["block_tag"]
+            },
+            artifacts={
+                "models" : exploration_steps.inputs.artifacts['models'],
+                "topology" :exploration_steps.inputs.artifacts['topology'],
+                "conf" : exploration_steps.inputs.artifacts['confs'],
+                "cv_file": exploration_steps.inputs.artifacts['cv_file']
+            },
+            key = step_keys["prep_exploration"]+"-{{item.order}}",
+            executor = prep_executor,
+            **prep_config
     )
     exploration_steps.add(prep_exploration)
 
-    run_exploration = Step(
+    run_merge = False
+    if run_executor is not None:
+        run_merge = run_executor.merge_sliced_step
+    if run_merge:
+        run_exploration = Step(
         'run-exploration',
         template=PythonOPTemplate(
             run_exploration_op,
@@ -171,7 +212,7 @@ def _exploration(
             retry_on_transient_error = retry_times,
             slices=Slices("{{item}}",
                 input_artifact=["task_path"],
-                output_artifact=["plm_out", "bias_fig","model_devi_fig", "trajectory", "md_log", "conf_out"]
+                output_artifact=["plm_out", "bias_fig","model_devi_fig", "dp_model_devi_fig", "dp_model_devi", "dp_selected_indices","dp_selected_confs","projected_fig","trajectory", "md_log", "conf_out"]
             ),
             **run_template_config,
         ),
@@ -190,8 +231,37 @@ def _exploration(
         key = step_keys["run_exploration"]+"-{{item}}",
         executor = run_executor,
         with_param=argo_range(argo_len(exploration_steps.inputs.parameters['task_names'])),
-        **run_config,
+        **run_config
     )
+    else:
+        run_exploration = Step(
+            'run-exploration',
+            template=PythonOPTemplate(
+                run_exploration_op,
+                python_packages = upload_python_package,
+                retry_on_transient_error = retry_times,
+                slices=Slices(sub_path = True,
+                    input_artifact=["task_path"],
+                    output_artifact=["plm_out", "bias_fig","model_devi_fig","dp_model_devi_fig", "dp_model_devi", "dp_selected_indices","dp_selected_confs","projected_fig","trajectory", "md_log", "conf_out"]
+                ),
+                **run_template_config,
+            ),
+            parameters={
+                "exploration_config" : exploration_steps.inputs.parameters["exploration_config"]
+            },
+            artifacts={
+                "task_path" : prep_exploration.outputs.artifacts["task_path"],
+                "forcefield": exploration_steps.inputs.artifacts['forcefield'],
+                "models" : exploration_steps.inputs.artifacts['models'],
+                "index_file": exploration_steps.inputs.artifacts['index_file'],
+                "dp_files": exploration_steps.inputs.artifacts['dp_files'],
+                "cv_file": exploration_steps.inputs.artifacts['cv_file'],
+                "inputfile": exploration_steps.inputs.artifacts['inputfile']
+            },
+            key = step_keys["run_exploration"]+"-{{item.order}}",
+            executor = run_executor,
+            **run_config,
+        )
     exploration_steps.add(run_exploration)
 
     exploration_steps.outputs.parameters["cv_dim"].value_from_parameter = prep_exploration.outputs.parameters["cv_dim"]
diff --git a/rid/superop/label.py b/rid/superop/label.py
index 48545ca8..a081a10e 100755
--- a/rid/superop/label.py
+++ b/rid/superop/label.py
@@ -30,6 +30,7 @@ def __init__(
         check_input_op: OP,
         prep_op: OP,
         run_op: OP,
+        stats_op: OP,
         prep_config: Dict,
         run_config: Dict,
         upload_python_package = None,
@@ -40,7 +41,7 @@ def __init__(
             "label_config": InputParameter(type=Dict),
             "cv_config": InputParameter(type=Dict),
             "tail": InputParameter(type=float, value=0.9),
-            "conf_tags": InputParameter(type=List),
+            "std_threshold": InputParameter(type=float, value=5.0),
             "block_tag" : InputParameter(type=str, value="")
         }        
         self._input_artifacts = {
@@ -52,7 +53,8 @@ def __init__(
             "at": InputArtifact(optional=True),
             "index_file": InputArtifact(optional=True),
             "dp_files": InputArtifact(optional=True),
-            "cv_file": InputArtifact(optional=True)
+            "cv_file": InputArtifact(optional=True),
+            "conf_tags": InputArtifact(optional=True)
         }
         self._output_parameters = {
         }
@@ -77,6 +79,7 @@ def __init__(
             "check_label_inputs": "{}-check-label-inputs".format(self.inputs.parameters["block_tag"]),
             "prep_label": "{}-prep-label".format(self.inputs.parameters["block_tag"]),
             "run_label": "{}-run-label".format(self.inputs.parameters["block_tag"]),
+            "label_stats": "{}-label-stats".format(self.inputs.parameters["block_tag"])
         }
 
         self = _label(
@@ -85,6 +88,7 @@ def __init__(
             check_input_op,
             prep_op,
             run_op,
+            stats_op,
             prep_config = prep_config,
             run_config = run_config,
             upload_python_package = upload_python_package,
@@ -118,6 +122,7 @@ def _label(
         check_label_input_op : OP,
         prep_label_op : OP,
         run_label_op : OP,
+        label_stats_op: OP,
         prep_config : Dict,
         run_config : Dict,
         upload_python_package : str = None,
@@ -140,10 +145,9 @@ def _label(
             retry_on_transient_error = retry_times,
             **prep_template_config,
         ),
-        parameters={
-            "conf_tags": label_steps.inputs.parameters['conf_tags'],  
-        },
+        parameters={},
         artifacts={
+            "conf_tags": label_steps.inputs.artifacts['conf_tags'],  
             "confs": label_steps.inputs.artifacts['confs'],
         },
         key = step_keys['check_label_inputs'],
@@ -152,7 +156,11 @@ def _label(
     )
     label_steps.add(check_label_inputs)
 
-    prep_label = Step(
+    prep_merge = False
+    if prep_executor is not None:
+        prep_merge = prep_executor.merge_sliced_step
+    if prep_merge:
+        prep_label = Step(
         'prep-label',
         template=PythonOPTemplate(
             prep_label_op,
@@ -179,11 +187,47 @@ def _label(
         executor = prep_executor,
         with_param=argo_range(argo_len(check_label_inputs.outputs.parameters['conf_tags'])),
         when = "%s > 0" % (check_label_inputs.outputs.parameters["if_continue"]),
-        **prep_config,
+        **prep_config
+    )
+    else:
+        prep_label = Step(
+        'prep-label',
+        template=PythonOPTemplate(
+            prep_label_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices("{{item}}",
+                group_size=10,
+                pool_size=1,
+                input_parameter=["task_name"],
+                input_artifact=["conf", "at"],
+                output_artifact=["task_path"]),
+            **prep_template_config,
+        ),
+        parameters={
+            "label_config": label_steps.inputs.parameters['label_config'],
+            "cv_config": label_steps.inputs.parameters['cv_config'],
+            "task_name": check_label_inputs.outputs.parameters['conf_tags']
+        },
+        artifacts={
+            "topology": label_steps.inputs.artifacts['topology'],
+            "conf": label_steps.inputs.artifacts['confs'],
+            "at": label_steps.inputs.artifacts['at'],
+            "cv_file": label_steps.inputs.artifacts['cv_file']
+        },
+        key = step_keys['prep_label']+"-{{item}}",
+        executor = prep_executor,
+        with_param=argo_range(argo_len(check_label_inputs.outputs.parameters['conf_tags'])),
+        when = "%s > 0" % (check_label_inputs.outputs.parameters["if_continue"]),
+        **prep_config
     )
     label_steps.add(prep_label)
 
-    run_label = Step(
+    run_merge = False
+    if run_executor is not None:
+        run_merge = run_executor.merge_sliced_step
+    if run_merge:
+        run_label = Step(
         'run-label',
         template=PythonOPTemplate(
             run_label_op,
@@ -192,7 +236,7 @@ def _label(
             slices=Slices("{{item}}",
                 input_parameter=["task_name"],
                 input_artifact=["task_path","at"],
-                output_artifact=["plm_out","cv_forces","mf_info","mf_fig","md_log"]),
+                output_artifact=["plm_out","cv_forces","mf_info","mf_fig","md_log","trajectory"]),
             **run_template_config,
         ),
         parameters={
@@ -214,11 +258,68 @@ def _label(
         executor = run_executor,
         with_param=argo_range(argo_len(check_label_inputs.outputs.parameters['conf_tags'])),
         continue_on_success_ratio = 0.75,
-        **run_config,
+        **run_config
+    )
+    else:
+        run_label = Step(
+        'run-label',
+        template=PythonOPTemplate(
+            run_label_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices("{{item}}",
+                group_size=10,
+                pool_size=1,
+                input_parameter=["task_name"],
+                input_artifact=["task_path","at"],
+                output_artifact=["plm_out","cv_forces","mf_info","mf_fig","md_log", "trajectory"]),
+            **run_template_config,
+        ),
+        parameters={
+            "label_config": label_steps.inputs.parameters["label_config"],
+            "cv_config": label_steps.inputs.parameters['cv_config'],
+            "task_name": check_label_inputs.outputs.parameters['conf_tags'],
+            "tail": label_steps.inputs.parameters['tail']
+        },
+        artifacts={
+            "forcefield": label_steps.inputs.artifacts['forcefield'],
+            "task_path": prep_label.outputs.artifacts["task_path"],
+            "index_file": label_steps.inputs.artifacts['index_file'],
+            "dp_files": label_steps.inputs.artifacts['dp_files'],
+            "cv_file": label_steps.inputs.artifacts['cv_file'],
+            "inputfile": label_steps.inputs.artifacts['inputfile'],
+            "at": label_steps.inputs.artifacts['at']
+        },
+        key = step_keys['run_label']+"-{{item}}",
+        executor = run_executor,
+        with_param=argo_range(argo_len(check_label_inputs.outputs.parameters['conf_tags'])),
+        continue_on_success_ratio = 0.75,
+        **run_config
     )
     label_steps.add(run_label)
+    
+    label_outputs_stats = Step(
+        'label-outputs-stats',
+        template=PythonOPTemplate(
+            label_stats_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            **prep_template_config,
+        ),
+        parameters={
+            "std_threshold": label_steps.inputs.parameters["std_threshold"]
+        },
+        artifacts={
+            "cv_forces": run_label.outputs.artifacts["cv_forces"],
+            "mf_info": run_label.outputs.artifacts["mf_info"]
+        },
+        key = step_keys["label_stats"],
+        executor = prep_executor,
+        **prep_config,
+    )
+    label_steps.add(label_outputs_stats)
 
-    label_steps.outputs.artifacts["cv_forces"]._from = run_label.outputs.artifacts["cv_forces"]
+    label_steps.outputs.artifacts["cv_forces"]._from = label_outputs_stats.outputs.artifacts["cv_forces"]
     label_steps.outputs.artifacts["md_log"]._from = run_label.outputs.artifacts["md_log"]
     
     return label_steps
diff --git a/rid/superop/mcmc.py b/rid/superop/mcmc.py
new file mode 100644
index 00000000..8e27351a
--- /dev/null
+++ b/rid/superop/mcmc.py
@@ -0,0 +1,200 @@
+from typing import Dict, List
+from copy import deepcopy
+from dflow import (
+    InputParameter,
+    OutputParameter,
+    Inputs,
+    InputArtifact,
+    Outputs,
+    OutputArtifact,
+    Step,
+    Steps,
+    argo_range,
+    argo_len,
+)
+from dflow.python import(
+    PythonOPTemplate,
+    OP,
+    Slices
+)
+from rid.utils import init_executor
+
+
+class MCMC(Steps):
+    
+    r"""" MCMC SuperOP.
+    This SuperOP combines MCMC_Run OP and MCMC_Plot OP.
+    """
+    def __init__(
+        self,
+        name: str,
+        mcmc_run_op: OP,
+        mcmc_plot_op: OP,
+        run_config: Dict,
+        plot_config: Dict,
+        upload_python_package = None,
+        retry_times = None
+    ):
+        self._input_parameters = {
+            "mcmc_config" : InputParameter(type=Dict),
+            "task_names" : InputParameter(type=List[str]),
+            "block_tag" : InputParameter(type=str, value="")
+        }        
+        self._input_artifacts = {
+            "models" : InputArtifact(),
+            "plm_out": InputArtifact()
+        }
+        self._output_parameters = {}
+        self._output_artifacts = {
+            "mcmc_fig": OutputArtifact()
+        }
+
+        super().__init__(        
+                name=name,
+                inputs=Inputs(
+                    parameters=self._input_parameters,
+                    artifacts=self._input_artifacts
+                ),
+                outputs=Outputs(
+                    parameters=self._output_parameters,
+                    artifacts=self._output_artifacts
+                ),
+            )
+
+        step_keys = {
+            "mcmc_run": "{}-mcmc-run".format(self.inputs.parameters["block_tag"]),
+            "mcmc_plot": "{}-mcmc-plot".format(self.inputs.parameters["block_tag"]),
+        }
+
+        self = _mcmc(
+            self, 
+            step_keys,
+            mcmc_run_op,
+            mcmc_plot_op,
+            run_config = run_config,
+            plot_config = plot_config,
+            upload_python_package = upload_python_package,
+            retry_times = retry_times
+        )            
+    
+    @property
+    def input_parameters(self):
+        return self._input_parameters
+
+    @property
+    def input_artifacts(self):
+        return self._input_artifacts
+
+    @property
+    def output_parameters(self):
+        return self._output_parameters
+
+    @property
+    def output_artifacts(self):
+        return self._output_artifacts
+
+    @property
+    def keys(self):
+        return self._keys
+
+
+def _mcmc(
+        mcmc_steps,
+        step_keys,
+        mcmc_run_op : OP,
+        mcmc_plot_op : OP,
+        run_config : Dict,
+        plot_config : Dict,
+        upload_python_package : str = None,
+        retry_times: int = None
+    ):
+    run_config = deepcopy(run_config)
+    plot_config = deepcopy(plot_config)
+    
+    run_template_config = run_config.pop('template_config')
+    plot_template_config = plot_config.pop('template_config')
+    run_executor = init_executor(run_config.pop('executor'))
+    plot_executor = init_executor(plot_config.pop('executor'))
+    
+    run_merge = False
+    if run_executor is not None:
+        run_merge = run_executor.merge_sliced_step
+    if run_merge:
+        mcmc_run = Step(
+        'mcmc-run',
+        template=PythonOPTemplate(
+            mcmc_run_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices("{{item}}",
+                input_parameter=["task_names"],
+                input_artifact=["models"],
+                output_artifact=["mcmc_1cv", "mcmc_2cv"]
+            ),
+            **run_template_config,
+        ),
+        parameters={
+            "mcmc_config" : mcmc_steps.inputs.parameters['mcmc_config'],
+            "task_names": mcmc_steps.inputs.parameters['task_names']
+        },
+        artifacts={
+            "models" : mcmc_steps.inputs.artifacts['models']
+        },
+        key = step_keys["mcmc_run"]+"-{{item}}",
+        with_param=argo_range(argo_len(mcmc_steps.inputs.parameters['task_names'])),
+        executor = run_executor,
+        **run_config
+    )
+    else:
+        mcmc_run = Step(
+            'mcmc-run',
+            template=PythonOPTemplate(
+                mcmc_run_op,
+                python_packages = upload_python_package,
+                retry_on_transient_error = retry_times,
+                slices=Slices(sub_path = True,
+                    input_parameter=["task_names"],
+                    input_artifact=["models"],
+                    output_artifact=["mcmc_1cv", "mcmc_2cv"]
+                ),
+                **run_template_config,
+            ),
+            parameters={
+                "mcmc_config" : mcmc_steps.inputs.parameters['mcmc_config'],
+                "task_names": mcmc_steps.inputs.parameters['task_names']
+            },
+            artifacts={
+                "models" : mcmc_steps.inputs.artifacts['models']
+            },
+            key = step_keys["mcmc_run"]+"-{{item.order}}",
+            executor = run_executor,
+            **run_config
+    )
+    mcmc_steps.add(mcmc_run)
+
+    mcmc_plot = Step(
+        'mcmc-plot',
+        template=PythonOPTemplate(
+            mcmc_plot_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            **plot_template_config,
+        ),
+        parameters={
+            "mcmc_config" : mcmc_steps.inputs.parameters["mcmc_config"]
+        },
+        artifacts={
+            "mcmc_1cv": mcmc_run.outputs.artifacts['mcmc_1cv'],
+            "mcmc_2cv": mcmc_run.outputs.artifacts['mcmc_2cv'],
+            "plm_out":  mcmc_steps.inputs.artifacts['plm_out']
+        },
+        key = step_keys["mcmc_plot"],
+        executor = plot_executor,
+        **run_config
+    )
+
+    mcmc_steps.add(mcmc_plot)
+
+    mcmc_steps.outputs.artifacts["mcmc_fig"]._from = mcmc_plot.outputs.artifacts["mcmc_fig"]
+    
+    return mcmc_steps
\ No newline at end of file
diff --git a/rid/superop/selector.py b/rid/superop/selector.py
index 8b819b0d..8a77c4b9 100755
--- a/rid/superop/selector.py
+++ b/rid/superop/selector.py
@@ -37,6 +37,7 @@ def __init__(
         retry_times = None
     ):
         self._input_parameters = {
+            "label_config": InputParameter(type=Dict),
             "trust_lvl_1" : InputParameter(type=List[float], value=2.0),
             "trust_lvl_2": InputParameter(type=List[float], value=3.0),
             "cluster_threshold": InputParameter(type=List[float], value=1.0),
@@ -48,6 +49,7 @@ def __init__(
             "dt": InputParameter(type=float, value=0.02),
             "output_freq": InputParameter(type=float, value=2500),
             "slice_mode": InputParameter(type=str, value="gmx"),
+            "type_map": InputParameter(type=List, value=[]),
             "if_make_threshold": InputParameter(type=bool, value=False),
             "task_names" : InputParameter(type=List[str]),
             "block_tag" : InputParameter(type=str, value="")
@@ -60,15 +62,15 @@ def __init__(
         }
         self._output_parameters = {
             "cluster_threshold": OutputParameter(type=List[int]),
-            "numb_cluster": OutputParameter(type=List[int]),
-            "selected_conf_tags": OutputParameter(type=List),
+            "numb_cluster": OutputParameter(type=List[int])
         }
         self._output_artifacts = {
             "cluster_selection_index": OutputArtifact(),
             "selected_confs": OutputArtifact(),
             "selected_cv_init": OutputArtifact(),
             "model_devi": OutputArtifact(),
-            "selected_indices": OutputArtifact()
+            "selected_indices": OutputArtifact(),
+            "selected_conf_tags": OutputArtifact()
         }
 
         super().__init__(        
@@ -138,7 +140,11 @@ def _select(
     prep_executor = init_executor(prep_config.pop('executor'))
     run_executor = init_executor(run_config.pop('executor'))
 
-    prep_select = Step(
+    prep_merge = False
+    if prep_executor is not None:
+        prep_merge = prep_executor.merge_sliced_step
+    if prep_merge:
+        prep_select = Step(
         'prep-select',
         template=PythonOPTemplate(
             prep_select_op,
@@ -170,11 +176,48 @@ def _select(
         key = step_keys["prep_select"]+"-{{item}}",
         executor = prep_executor,
         with_param=argo_range(argo_len(select_steps.inputs.parameters['task_names'])),
+        **prep_config
+    )  
+    else:
+        prep_select = Step(
+        'prep-select',
+        template=PythonOPTemplate(
+            prep_select_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices(sub_path = True,
+                input_parameter=["cluster_threshold", "task_name"],
+                input_artifact=["plm_out"],
+                output_artifact=["cluster_fig","cluster_selection_index", "cluster_selection_data"],
+                output_parameter=["cluster_threshold", "numb_cluster"]
+                ),
+                
+            **prep_template_config,
+        ),
+        parameters={
+            "cluster_threshold": select_steps.inputs.parameters['cluster_threshold'],
+            "angular_mask": select_steps.inputs.parameters['angular_mask'],
+            "weights": select_steps.inputs.parameters['weights'],
+            "numb_cluster_upper": select_steps.inputs.parameters['numb_cluster_upper'],
+            "numb_cluster_lower": select_steps.inputs.parameters['numb_cluster_lower'],
+            "max_selection": select_steps.inputs.parameters['max_selection'],
+            "if_make_threshold": select_steps.inputs.parameters['if_make_threshold'],
+            "task_name": select_steps.inputs.parameters['task_names']
+        },
+        artifacts={
+            "plm_out": select_steps.inputs.artifacts['plm_out']
+        },
+        key = step_keys["prep_select"]+"-{{item.order}}",
+        executor = prep_executor,
         **prep_config,
     )
     select_steps.add(prep_select)
 
-    run_select = Step(
+    run_merge = False
+    if run_executor is not None:
+        run_merge = run_executor.merge_sliced_step
+    if run_merge:
+        run_select = Step(
         'run-select',
         template=PythonOPTemplate(
             run_select_op,
@@ -184,17 +227,18 @@ def _select(
                 "int({{item}})",
                 input_parameter=["task_name", "trust_lvl_1", "trust_lvl_2"],
                 input_artifact=["cluster_selection_index", "cluster_selection_data", "xtc_traj", "topology"],
-                output_artifact=["selected_confs", "selected_cv_init", "model_devi", "selected_indices"],
-                output_parameter=["selected_conf_tags"]
+                output_artifact=["selected_confs", "selected_cv_init", "model_devi", "selected_indices","selected_conf_tags"]
             ),
             **run_template_config,
         ),
         parameters={
+            "label_config":  select_steps.inputs.parameters["label_config"],
             "trust_lvl_1": select_steps.inputs.parameters["trust_lvl_1"],
             "trust_lvl_2": select_steps.inputs.parameters["trust_lvl_2"],
             "dt": select_steps.inputs.parameters["dt"],
             "output_freq": select_steps.inputs.parameters["output_freq"],
             "slice_mode": select_steps.inputs.parameters["slice_mode"],
+            "type_map": select_steps.inputs.parameters["type_map"],
             "task_name": select_steps.inputs.parameters['task_names']
         },
         artifacts={
@@ -207,14 +251,49 @@ def _select(
         key = step_keys["run_select"]+"-{{item}}",
         executor = run_executor,
         with_param=argo_range(argo_len(select_steps.inputs.parameters["task_names"])),
+        **run_config
+    )
+    else:
+        run_select = Step(
+        'run-select',
+        template=PythonOPTemplate(
+            run_select_op,
+            python_packages = upload_python_package,
+            retry_on_transient_error = retry_times,
+            slices=Slices(sub_path = True,
+                input_parameter=["task_name", "trust_lvl_1", "trust_lvl_2"],
+                input_artifact=["cluster_selection_index", "cluster_selection_data", "xtc_traj", "topology"],
+                output_artifact=["selected_confs", "selected_cv_init", "model_devi", "selected_indices","selected_conf_tags"]
+            ),
+            **run_template_config,
+        ),
+        parameters={
+            "label_config":  select_steps.inputs.parameters["label_config"],
+            "trust_lvl_1": select_steps.inputs.parameters["trust_lvl_1"],
+            "trust_lvl_2": select_steps.inputs.parameters["trust_lvl_2"],
+            "dt": select_steps.inputs.parameters["dt"],
+            "output_freq": select_steps.inputs.parameters["output_freq"],
+            "slice_mode": select_steps.inputs.parameters["slice_mode"],
+            "type_map": select_steps.inputs.parameters["type_map"],
+            "task_name": select_steps.inputs.parameters['task_names']
+        },
+        artifacts={
+            "cluster_selection_index": prep_select.outputs.artifacts["cluster_selection_index"],
+            "cluster_selection_data": prep_select.outputs.artifacts["cluster_selection_data"],
+            "models": select_steps.inputs.artifacts["models"],
+            "xtc_traj": select_steps.inputs.artifacts["xtc_traj"],
+            "topology": select_steps.inputs.artifacts["topology"]
+        },
+        key = step_keys["run_select"]+"-{{item.order}}",
+        executor = run_executor,
         **run_config,
     )
     select_steps.add(run_select)
 
     select_steps.outputs.parameters["cluster_threshold"].value_from_parameter = prep_select.outputs.parameters["cluster_threshold"]
     select_steps.outputs.parameters["numb_cluster"].value_from_parameter = prep_select.outputs.parameters["numb_cluster"]
-    select_steps.outputs.parameters["selected_conf_tags"].value_from_parameter = run_select.outputs.parameters["selected_conf_tags"]
 
+    select_steps.outputs.artifacts["selected_conf_tags"]._from = run_select.outputs.artifacts["selected_conf_tags"]
     select_steps.outputs.artifacts["cluster_selection_index"]._from = prep_select.outputs.artifacts["cluster_selection_index"]
     select_steps.outputs.artifacts["selected_confs"]._from = run_select.outputs.artifacts["selected_confs"]
     select_steps.outputs.artifacts["selected_cv_init"]._from = run_select.outputs.artifacts["selected_cv_init"]
diff --git a/rid/template/machine_bohrium.json b/rid/template/machine_bohrium.json
index 888069de..95e4e053 100755
--- a/rid/template/machine_bohrium.json
+++ b/rid/template/machine_bohrium.json
@@ -2,13 +2,13 @@
     "resources": {
         "bohrium1": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -32,13 +32,13 @@
         },
         "bohrium2": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -62,13 +62,13 @@
         },
         "bohrium3": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -92,13 +92,13 @@
         },
         "bohrium4": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -122,13 +122,13 @@
         },
         "bohrium5": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -161,6 +161,7 @@
         "run_select_config": "bohrium4",
         "prep_data_config": "bohrium3",
         "run_train_config": "bohrium5",
+        "model_devi_config": "bohrium5",
         "workflow_steps_config": "bohrium3"
     }
 }
\ No newline at end of file
diff --git a/rid/template/machine_bohrium_k8s.json b/rid/template/machine_bohrium_k8s.json
index b4378035..8fc03b57 100755
--- a/rid/template/machine_bohrium_k8s.json
+++ b/rid/template/machine_bohrium_k8s.json
@@ -4,18 +4,18 @@
             "template_config" : {
                 "image": "registry.dp.tech/public/pkufjhdocker/rid-tf-cpu:stable", 
                 "image_pull_policy" : "IfNotPresent",
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             }
         },
         "bohrium1": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -39,13 +39,13 @@
         },
         "bohrium2": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -69,13 +69,13 @@
         },
         "bohrium3": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -99,13 +99,13 @@
         },
         "bohrium4": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
             "machine_dict":{
-                "batch_type": "DpCloudServer",
-                "context_type": "DpCloudServerContext",
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
                 "local_root" : "./",
                 "remote_profile":{
                     "email": "",
@@ -138,6 +138,7 @@
         "run_select_config": "bohrium3",
         "prep_data_config": "local_k8s",
         "run_train_config": "bohrium4",
+        "model_devi_config": "bohrium4",
         "workflow_steps_config": "local_k8s"
     }
 }
\ No newline at end of file
diff --git a/rid/template/machine_bohrium_k8s_mcmc.json b/rid/template/machine_bohrium_k8s_mcmc.json
new file mode 100644
index 00000000..2c4b90a0
--- /dev/null
+++ b/rid/template/machine_bohrium_k8s_mcmc.json
@@ -0,0 +1,79 @@
+{
+    "resources": {
+        "local_k8s": {
+            "template_config" : {
+                "image": "registry.dp.tech/public/pkufjhdocker/rid-tf-cpu:stable", 
+                "image_pull_policy" : "IfNotPresent",
+                "requests" : {"ephemeral-storage": "4Gi"}
+            }
+        },
+        "bohrium1": {
+            "template_config":{
+                "requests" : {"ephemeral-storage": "4Gi"}
+            },
+            "executor":{
+            "image": "dptechnology/dpdispatcher:latest",
+            "machine_dict":{
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
+                "local_root" : "./",
+                "remote_profile":{
+                    "email": "",
+                    "password": "",
+                    "program_id": "",
+                    "input_data":{
+                        "api_version":2,
+                        "job_type": "container",
+                        "job_group_id": "",
+                        "log_file": "tmp_log",
+                        "job_name": "train",
+                        "scass_type":"c12_m92_1 * NVIDIA V100",
+                        "platform": "ali",
+                        "image_name":"registry.dp.tech/public/pkufjhdocker/rid-tf-gpu:stable"
+                        }
+                }
+            },
+            "resources_dict":{
+                "batch_type:": "Bohrium"
+            }
+        }
+        },
+        "bohrium2": {
+            "template_config":{
+                "requests" : {"ephemeral-storage": "4Gi"}
+            },
+            "executor":{
+            "image": "dptechnology/dpdispatcher:latest",
+            "machine_dict":{
+                "batch_type": "Bohrium",
+                "context_type": "Bohrium",
+                "local_root" : "./",
+                "remote_profile":{
+                    "email": "",
+                    "password": "",
+                    "program_id": "",
+                    "input_data":{
+                        "api_version":2,
+                        "job_type": "container",
+                        "job_group_id": "",
+                        "log_file": "tmp_log",
+                        "job_name": "train",
+                        "scass_type":"c12_m92_1 * NVIDIA V100",
+                        "platform": "ali",
+                        "image_name":"registry.dp.tech/public/pkufjhdocker/rid-tf-cpu:stable"
+                        }
+                }
+            },
+            "resources_dict":{
+                "batch_type:": "Bohrium"
+            }
+        }
+        }
+    },
+
+    "tasks": {
+        "mcmc_run_config": "bohrium1",
+        "mcmc_plot_config": "bohrium2",
+        "workflow_steps_config": "local_k8s"
+    }
+}
\ No newline at end of file
diff --git a/rid/template/machine_local_k8s.json b/rid/template/machine_local_k8s.json
index e8dbefe2..43b628d1 100644
--- a/rid/template/machine_local_k8s.json
+++ b/rid/template/machine_local_k8s.json
@@ -46,6 +46,7 @@
         "run_select_config": "local_k8s_4",
         "prep_data_config": "local_k8s_3",
         "run_train_config": "local_k8s_5",
+        "model_devi_config": "local_k8s_5",
         "workflow_steps_config": "local_k8s_3"
     }
 }
\ No newline at end of file
diff --git a/rid/template/machine_slurm_k8s.json b/rid/template/machine_slurm_k8s.json
index b23d36e7..b55efefa 100644
--- a/rid/template/machine_slurm_k8s.json
+++ b/rid/template/machine_slurm_k8s.json
@@ -4,12 +4,12 @@
             "template_config" : {
                 "image": "registry.dp.tech/public/pkufjhdocker/rid-tf-cpu:stable", 
                 "image_pull_policy" : "IfNotPresent",
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             }
         },
         "remote_slurm": {
             "template_config":{
-                "requests" : {"ephemeral-storage": "1Gi"}
+                "requests" : {"ephemeral-storage": "4Gi"}
             },
             "executor":{
             "image": "dptechnology/dpdispatcher:latest",
@@ -51,6 +51,7 @@
         "run_select_config": "remote_slurm",
         "prep_data_config": "local_k8s",
         "run_train_config": "remote_slurm",
+        "model_devi_config": "remote_slurm",
         "workflow_steps_config": "local_k8s"
     }
 }
\ No newline at end of file
diff --git a/rid/template/machine_slurm_k8s_mcmc.json b/rid/template/machine_slurm_k8s_mcmc.json
new file mode 100644
index 00000000..8be5069c
--- /dev/null
+++ b/rid/template/machine_slurm_k8s_mcmc.json
@@ -0,0 +1,50 @@
+{
+    "resources": {
+        "local_k8s": {
+            "template_config" : {
+                "image": "registry.dp.tech/public/pkufjhdocker/rid-tf-cpu:stable", 
+                "image_pull_policy" : "IfNotPresent",
+                "requests" : {"ephemeral-storage": "4Gi"}
+            }
+        },
+        "remote_slurm": {
+            "template_config":{
+                "requests" : {"ephemeral-storage": "4Gi"}
+            },
+            "executor":{
+            "image": "dptechnology/dpdispatcher:latest",
+            "machine_dict":{
+                "batch_type": "Slurm",
+                "context_type": "SSHContext",
+                "local_root" : "./",
+                "remote_root": "",
+                "remote_profile":{
+                    "hostname": "",
+                    "username": "",
+                    "password": "",
+                    "port": "",
+                    "timeout": 20
+                }
+            },
+            "resources_dict":{
+                "number_node": 1,
+                "cpu_per_node": 8,
+                "gpu_per_node": 1,
+                "queue_name": "",
+                "group_size": 1,
+                "custom_flags": [
+                    "#SBATCH --time=120:00:00",
+                    "#SBATCH --exclude=3090-[000,007]"
+                   ],
+                "source_list": ["/path/to/rid.env"]
+            }
+        }
+        }
+    },
+
+    "tasks": {
+        "mcmc_run_config": "remote_slurm",
+        "mcmc_plot_config": "remote_slurm",
+        "workflow_steps_config": "local_k8s"
+    }
+}
\ No newline at end of file
diff --git a/rid/template/machine_slurm_local.json b/rid/template/machine_slurm_local.json
index fc29c60e..9c22bd87 100644
--- a/rid/template/machine_slurm_local.json
+++ b/rid/template/machine_slurm_local.json
@@ -57,6 +57,7 @@
         "run_select_config": "local_machine",
         "prep_data_config": "local_machine",
         "run_train_config": "local_machine",
+        "model_devi_config": "local_machine",
         "workflow_steps_config": "local_machine_cpu"
     }
 }
\ No newline at end of file
diff --git a/rid/template/machine_slurm_local_mcmc.json b/rid/template/machine_slurm_local_mcmc.json
new file mode 100644
index 00000000..8e5c9f6a
--- /dev/null
+++ b/rid/template/machine_slurm_local_mcmc.json
@@ -0,0 +1,56 @@
+{
+    "resources": {
+        "local_machine": {
+            "executor":{
+            "image": "dptechnology/dpdispatcher:latest",
+            "merge_sliced_step": "True",
+            "machine_dict":{
+                "batch_type": "Slurm",
+                "context_type": "Local",
+                "local_root" : "./",
+                "remote_root": ""
+            },
+            "resources_dict":{
+                "number_node": 1,
+                "cpu_per_node": 8,
+                "gpu_per_node": 1,
+                "queue_name": "GPU",
+                "group_size": 1,
+                "custom_flags": [
+                    "#SBATCH --time=120:00:00"
+                   ],
+                "source_list": ["/path/to/rid.env"]
+            }
+        }
+        },
+        "local_machine_cpu": {
+            "executor":{
+            "image": "dptechnology/dpdispatcher:latest",
+            "merge_sliced_step": "True",
+            "machine_dict":{
+                "batch_type": "Slurm",
+                "context_type": "Local",
+                "local_root" : "./",
+                "remote_root": ""
+            },
+            "resources_dict":{
+                "number_node": 1,
+                "cpu_per_node": 8,
+                "gpu_per_node": 0,
+                "queue_name": "C032M0128G",
+                "group_size": 1,
+                "custom_flags": [
+                    "#SBATCH --time=120:00:00"
+                   ],
+                "source_list": ["/path/to/rid.env"]
+            }
+        }
+        }
+    },
+
+    "tasks": {
+        "mcmc_run_config": "remote_slurm",
+        "mcmc_plot_config": "remote_slurm",
+        "workflow_steps_config": "local_k8s"
+    }
+}
\ No newline at end of file
diff --git a/rid/template/rid_gmx_dih.json b/rid/template/rid_gmx_dih.json
index 7f5524b0..88c49468 100755
--- a/rid/template/rid_gmx_dih.json
+++ b/rid/template/rid_gmx_dih.json
@@ -1,7 +1,7 @@
 {
     "name": "test",
-    "numb_walkers": 2,
-    "numb_iters": 2,
+    "numb_walkers": 4,
+    "numb_iters": 5,
     "trust_lvl_1": 2,
     "trust_lvl_2": 3,
     "init_models": [],
@@ -18,7 +18,7 @@
         "nsteps": 50000,
         "type": "gmx",
         "temperature": 300,
-        "output_freq": 500,
+        "output_freq": 50,
         "ref-t": "300 300",
         "verlet-buffer-tolerance":"-1",
         "rlist": 1,
@@ -50,7 +50,7 @@
         "temperature":300,
         "method": "restrained",
         "type": "gmx",
-        "output_freq": 100,
+        "output_freq": 50,
         "ref-t": "300 300",
         "rlist": 1,
         "verlet-buffer-tolerance":"-1",
@@ -66,19 +66,20 @@
         "ntmpi": 1,
         "nt": 8,
         "max_warning": 2,
-        "kappas": [ 500, 500 ]
+        "kappas": [ 500, 500 ],
+        "std_threshold": 1.5
     },
 
     "Train": {
         "numb_models": 4,
-        "neurons": [ 200, 200, 200, 200 ],
+        "neurons": [40, 40, 40, 40],
         "resnet": true,
         "batch_size": 128,
         "epoches": 8000,
         "init_lr": 0.0008,
         "decay_steps": 120,
         "decay_rate": 0.96,
-        "drop_out_rate": 0,
+        "drop_out_rate": 0.1,
         "numb_threads": 8,
         "use_mix": false,
         "restart": false
diff --git a/rid/template/rid_gmx_dis.json b/rid/template/rid_gmx_dis.json
index b2a747fe..7b1311e0 100755
--- a/rid/template/rid_gmx_dis.json
+++ b/rid/template/rid_gmx_dis.json
@@ -66,7 +66,8 @@
         "output_mode": "both",
         "ntmpi": 1,
         "nt": 8,
-        "max_warning": 2
+        "max_warning": 2,
+        "std_threshold": 5.0
     },
 
     "Train": {
@@ -78,7 +79,7 @@
         "init_lr": 0.0008,
         "decay_steps": 120,
         "decay_rate": 0.96,
-        "drop_out_rate": 0,
+        "drop_out_rate": 0.1,
         "numb_threads": 8,
         "use_mix": false,
         "restart": false
diff --git a/rid/template/rid_gmx_dp.json b/rid/template/rid_gmx_dp.json
index a43acc9f..0cacdae9 100755
--- a/rid/template/rid_gmx_dp.json
+++ b/rid/template/rid_gmx_dp.json
@@ -81,7 +81,8 @@
         "nstenergy": 100,
         "nt": 8,
         "max_warning": 3,
-        "kappas": [500, 500]
+        "kappas": [500, 500],
+        "std_threshold": 5.0
     },
 
     "Train": {
diff --git a/rid/template/rid_lmp_dp.json b/rid/template/rid_lmp_dp.json
index 3f9cd2dd..b756c345 100755
--- a/rid/template/rid_lmp_dp.json
+++ b/rid/template/rid_lmp_dp.json
@@ -43,7 +43,8 @@
         "kappas": [ 500, 500, 500, 500],
         "ntmpi": 1,
         "nt": 8,
-        "max_warning": 2
+        "max_warning": 2,
+        "std_threshold": 1.5
     },
 
     "Train": {
diff --git a/rid/template/rid_lmp_gmx_dp.json b/rid/template/rid_lmp_gmx_dp.json
new file mode 100755
index 00000000..6f679568
--- /dev/null
+++ b/rid/template/rid_lmp_gmx_dp.json
@@ -0,0 +1,101 @@
+{
+    "name": "rid-mixed-gmx-lmp",
+    "numb_walkers": 2,
+    "numb_iters": 2,
+    "trust_lvl_1": 2,
+    "trust_lvl_2": 4,
+    "init_models": [],
+    
+    "CV": {
+        "mode": "custom",
+        "selected_atomid":[[161,165],[124,156],[124,161],[156,165],[107,161],[139,161],[153,161],[121,162],[121,163],[121,164],[138,161],[154,161],[153,165]],
+        "angular_mask": [0,0,0,0,0,0,0,0,0,0,0,0,0],
+        "weights": [1,1,1,1,1,1,1,1,1,1,1,1,1],
+        "units": "A",
+        "cv_file": ["colvar", "plmpath.pdb"]
+    },
+
+    "SelectorConfig": {
+        "numb_cluster_lower": 6,
+        "numb_cluster_upper": 10,
+        "cluster_threshold": 2.0,
+        "max_selection": 20,
+        "numb_cluster_threshold": 5,
+        "slice_mode": "dpdata",
+        "type_map": ["Si", "O", "H", "Hf", "Cl"]
+    },
+
+    "ExploreMDConfig": {
+        "type": "lmp",
+        "temperature": 372,
+        "dp_model": ["1a0.pb","1a1.pb","1a2.pb","1a3.pb"],
+        "dp_model_devi_out": "model_devi.out",
+        "dp_model_devi_f_lo": 0.14,
+        "dp_model_devi_f_hi": 0.26,
+        "dp_model_max_selection": 20,
+        "proj_info": {
+            "proj_mode": "path",
+            "proj_cv_index": [0,1,2,3],
+            "path_list": [[2.38722,0.99005,3.99886,2.08580],[2.64067,1.20917,3.73831,1.86004],[2.90639,1.50430,3.47935,1.65185],
+            [3.18130,1.83919,3.22240,1.46871],[3.46322,2.19574,2.96795,1.32108],[3.75057,2.56491,2.71672,1.22191],
+            [4.04218,2.94198,2.46968,1.18344],[4.33720,3.32424,2.22824,1.21147],[4.63498,3.71010,1.99443,1.30170]],
+            "path_lm": 0.25
+        },
+        "inputfile":"input_explore.lammps",
+        "dt": 0.001,
+        "nsteps": 10000,
+        "output_freq": 100,
+        "ntmpi": 1,
+        "nt": 12,
+        "max_warning": 2
+    },
+
+    "LabelMDConfig": {
+        "type": "gmx",
+        "method": "constrained",
+        "temperature": 372,
+        "dp_model": ["1a0.pb"],
+        "nsteps": 10000,
+        "ref-t": "372 0",
+        "gen-vel": "yes",
+        "gen-temp": "372",
+        "tc-grps": "top bottom",
+        "dt": 0.001,
+        "tau-t": "0.05 0.05",
+        "constraints": "none", 
+        "comm-grps": "top",
+        "output_freq": 100,
+        "output_mode": "both",
+        "freezegrps": "bottom",
+        "freezedim": "Y Y Y",
+        "pcoupl": "no",
+        "ntmpi": 1,
+        "coulombtype": "cutoff",
+        "rlist": 1,
+        "rvdw": 0.7,
+        "rvdw-switch": 0,
+        "rcoulomb": 0.7,
+        "rcoulomb-switch": 0,
+        "nstenergy": 1000,
+        "nt": 12,
+        "max_warning": 3,
+        "std_threshold": 5.0
+    },
+
+    "Train": {
+        "numb_models": 4,
+        "neurons": [120, 120, 120, 120],
+        "resnet": true,
+        "batch_size": 4096,
+        "epoches": 30000,
+        "init_lr": 0.001,
+        "decay_steps": 200,
+        "decay_rate": 0.96,
+        "train_thread": 8,
+        "drop_out_rate": 0.1,
+        "numb_threads": 8,
+        "use_mix": false,
+        "restart": false,
+        "decay_steps_inner": 200
+    }
+}
\ No newline at end of file
diff --git a/rid/template/rid_mcmc_cv_dih.json b/rid/template/rid_mcmc_cv_dih.json
new file mode 100755
index 00000000..e2cbc6a1
--- /dev/null
+++ b/rid/template/rid_mcmc_cv_dih.json
@@ -0,0 +1,18 @@
+{
+        "name": "rid-mcmc",
+        "init_models": ["model_000.pb","model_001.pb","model_002.pb","model_003.pb"],
+        
+        "MCMC_Config": {
+            "cv_dimension": 2,
+            "numb_steps": 10000,
+            "numb_walkers": 500,
+            "temperature": 300,
+            "proj_info": {
+                "proj_mode": "cv",
+                "proj_cv_index": [0,1]
+            },
+            "cv_type": "dih",
+            "bins": 101
+        }
+    
+}
\ No newline at end of file
diff --git a/rid/template/rid_mcmc_cv_dis.json b/rid/template/rid_mcmc_cv_dis.json
new file mode 100755
index 00000000..5f0a707d
--- /dev/null
+++ b/rid/template/rid_mcmc_cv_dis.json
@@ -0,0 +1,20 @@
+{
+        "name": "rid-mcmc",
+        "init_models": ["model_000.pb","model_001.pb","model_002.pb","model_003.pb"],
+        
+        "MCMC_Config": {
+            "cv_dimension": 20,
+            "numb_steps": 100000,
+            "numb_walkers": 2000,
+            "temperature": 300,
+            "cv_upper_bound": [4.5,5.0,6.0,4.2,4.0,4.0,3.5,4.0,4.0,3.5,3.2,3.0,3.5,3.5,4.0,4.5,4.5,3.75,3.8,4.0],
+            "cv_lower_bound": [1.0,2.0,3.5,2.8,2.5,2.0,1.2,1.5,1.5,1.0,1.8,0.5,1.0,1.0,1.5,1.5,1.5,1.75,2.2,2.6],
+            "proj_info": {
+                "proj_mode": "cv",
+                "proj_cv_index": [0,1]
+            },
+            "cv_type": "dis",
+            "bins": 101
+        }
+    
+}
\ No newline at end of file
diff --git a/rid/template/rid_mcmc_path.json b/rid/template/rid_mcmc_path.json
new file mode 100755
index 00000000..451ea58f
--- /dev/null
+++ b/rid/template/rid_mcmc_path.json
@@ -0,0 +1,24 @@
+{
+        "name": "rid-mcmc",
+        "init_models": ["model_000.pb","model_001.pb","model_002.pb","model_003.pb"],
+        
+        "MCMC_Config": {
+            "cv_dimension": 13,
+            "numb_steps": 100000,
+            "numb_walkers": 2000,
+            "temperature": 300,
+            "cv_upper_bound": [5.33507,4.43402,4.59772,3.40287,4.67268,5.12813,3.28143,3.74967,3.58144,5.34143,7.13014,4.58024,7.41678],
+            "cv_lower_bound": [2.38722,0.99005,3.99886,2.08580,3.81799,3.93331,2.768175,3.145618,3.161886,4.661101,5.286805,4.076754,2.893047],
+            "proj_info": {
+                "proj_mode": "path",
+                "proj_cv_index": [0,1,2,3],
+                "path_list": [[2.38722,0.99005,3.99886,2.08580],[2.64067,1.20917,3.73831,1.86004],[2.90639,1.50430,3.47935,1.65185],
+                [3.18130,1.83919,3.22240,1.46871],[3.46322,2.19574,2.96795,1.32108],[3.75057,2.56491,2.71672,1.22191],
+                [4.04218,2.94198,2.46968,1.18344],[4.33720,3.32424,2.22824,1.21147],[4.63498,3.71010,1.99443,1.30170]],
+                "path_lm": 0.25
+            },
+            "cv_type": "dis",
+            "bins": 101
+        }
+    
+}
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 580eaa3a..da683539 100644
--- a/setup.py
+++ b/setup.py
@@ -14,10 +14,11 @@
 readme = readme_file.read_text(encoding="utf-8")
 
 setuptools.setup(
-    name='rid',
+    name='rid-kit',
     author="Yanze Wang, Jiahao Fan",
+    author_email="yanze039@mit.edu,jiahaofan@pku.edu.cn",
     use_scm_version={'write_to': 'rid/_version.py'},
-    author_email="yanze039@mit.edu;jiahaofan@pku.edu.cn",
+    version='{{VERSION_PLACEHOLDER}}',
     description="RiD package for enhanced sampling",
     setup_requires=['setuptools_scm'],
     long_description=readme,
@@ -36,6 +37,7 @@
         "rid/nn",
         "rid/op",
         "rid/select",
+        "rid/mcmc",
         "rid/superop",
         "rid/task",
         "rid/utils",


Test patch:
diff --git a/tests/data/conf.json b/tests/data/conf.json
new file mode 100644
index 00000000..03eeec12
--- /dev/null
+++ b/tests/data/conf.json
@@ -0,0 +1,1 @@
+{"conf.gro":"000"}
\ No newline at end of file
diff --git a/tests/data/model_000.pb b/tests/data/model_000.pb
deleted file mode 100644
index eedf9dc3..00000000
Binary files a/tests/data/model_000.pb and /dev/null differ
diff --git a/tests/data/model_001.pb b/tests/data/model_001.pb
deleted file mode 100644
index a06c0362..00000000
Binary files a/tests/data/model_001.pb and /dev/null differ
diff --git a/tests/data/model_002.pb b/tests/data/model_002.pb
deleted file mode 100644
index 368902ff..00000000
Binary files a/tests/data/model_002.pb and /dev/null differ
diff --git a/tests/data/models/model_000.pb b/tests/data/models/model_000.pb
new file mode 100644
index 00000000..1ed33b56
Binary files /dev/null and b/tests/data/models/model_000.pb differ
diff --git a/tests/data/models/model_001.pb b/tests/data/models/model_001.pb
new file mode 100644
index 00000000..f44b3b5d
Binary files /dev/null and b/tests/data/models/model_001.pb differ
diff --git a/tests/data/models/model_002.pb b/tests/data/models/model_002.pb
new file mode 100644
index 00000000..f439d618
Binary files /dev/null and b/tests/data/models/model_002.pb differ
diff --git a/tests/data/models/model_003.pb b/tests/data/models/model_003.pb
new file mode 100644
index 00000000..ca409ddf
Binary files /dev/null and b/tests/data/models/model_003.pb differ
diff --git a/tests/data/rid.json b/tests/data/rid.json
index 10279188..317e307d 100755
--- a/tests/data/rid.json
+++ b/tests/data/rid.json
@@ -25,6 +25,7 @@
 
     "ExploreMDConfig": {
         "nsteps": 50,
+        "type": "gmx",
         "output_freq": 1,
         "temperature": 300,
         "dt": 0.002,
@@ -36,6 +37,7 @@
 
     "LabelMDConfig": {
         "nsteps": 50,
+        "type": "gmx",
         "output_freq": 1,
         "temperature": 300,
         "dt": 0.002,
@@ -43,7 +45,8 @@
         "ntmpi": 1,
         "nt": 8,
         "max_warning": 0,
-        "kappas": [ 500, 500 ]
+        "kappas": [ 500, 500 ],
+        "std_threshold": 2.0
     },
 
     "Train": {
diff --git a/tests/mocked_ops.py b/tests/mocked_ops.py
index 686edb0c..17fd4959 100755
--- a/tests/mocked_ops.py
+++ b/tests/mocked_ops.py
@@ -23,12 +23,13 @@
 from rid.op.prep_label import CheckLabelInputs
 from rid.op.run_exploration import RunExplore
 from rid.op.run_label import RunLabel
+from rid.op.label_stats import LabelStats
 from rid.op.prep_select import PrepSelect
 from rid.op.run_select import RunSelect
 from rid.op.run_train import TrainModel
 from rid.constants import (
     tf_model_name,
-    init_conf_name,
+    init_conf_gmx_name,
     data_new,
     data_raw,
     gmx_conf_name,
@@ -47,7 +48,8 @@
     cv_init_label,
     model_devi_precision,
     train_log,
-    train_fig
+    train_fig,
+    mf_std_fig
 )
 
 upload_packages.append(__file__)
@@ -59,9 +61,13 @@ def clear_dir(data=None):
 
 def clear_files(files=None):
     if files:
-        for file in files:
-            if os.path.exists(file):
-                os.remove(file)
+        if isinstance(files,list):
+            for file in files:
+                if os.path.exists(file):
+                    os.remove(file)
+        elif isinstance(files,Path):
+            if os.path.exists(files):
+                    os.remove(files)
 
 def make_mocked_init_models(numb_models):
     tmp_models = []
@@ -85,7 +91,7 @@ def make_mocked_init_data(data_name,data_num):
 def make_mocked_init_confs(numb_confs):
     tmp_confs = []
     for ii in range(numb_confs):
-        ff = Path(init_conf_name.format(idx=ii))
+        ff = Path(init_conf_gmx_name.format(idx=ii))
         ff.write_text(f'This is init conf {ii}')
         tmp_confs.append(ff)
     return tmp_confs
@@ -95,6 +101,12 @@ def make_mocked_file(file_name):
     ff.write_text(f'This is mocked %s'%file_name)
     return ff
 
+def make_mocked_json(filename,fconts):
+    ff = Path(filename)
+    with open(ff,"w") as f:
+        json.dump(fconts,f)
+    return ff
+
 class MockedCollectData(CollectData):
     @OP.exec_sign_check
     def execute(
@@ -164,6 +176,7 @@ def execute(
             self,
             ip : OPIO,
     ) -> OPIO:
+        dp_conf_list = []
         assert(ip["task_path"].is_dir())
         print("start run explore !!")
         with set_directory(ip["task_path"]):
@@ -183,6 +196,7 @@ def execute(
                 f.write("gmx trajectory file")
         
         op = OPIO({
+                "dp_selected_confs": dp_conf_list,
                 "plm_out": ip["task_path"].joinpath(plumed_output_name),
                 "md_log": ip["task_path"].joinpath(gmx_mdrun_log),
                 "trajectory": ip["task_path"].joinpath(gmx_xtc_name),
@@ -204,12 +218,12 @@ def execute(
 
             tags = {}
             for tag in ip["conf_tags"]:
-                if isinstance(tag, Dict):
-                    tags.update(tag)
-                elif isinstance(tag, str):
-                    tags.update(json.loads(tag))
+                if isinstance(tag,Path):
+                    with open(tag,"r") as f:
+                        tags.update(json.load(f))
                 else:
                     raise RuntimeError("Unkown Error.")
+                
             conf_tags = []
             for conf in ip["confs"]:
                 conf_tags.append(str(tags[conf.name]))
@@ -260,13 +274,36 @@ def execute(
                 f.write("gmx md log file")
             with open(cv_force_out,"w") as f:
                 f.write('1.000000e+00 2.000000e+00\n')
+            with open("mf_info.out", "w") as f:
+                f.write("mean force information file")
+            with open(gmx_xtc_name, "w") as f:
+                f.write("gmx trajectory")
         
         op = OPIO({
                 "plm_out": ip["task_path"].joinpath(plumed_output_name),
+                "trajectory": ip["task_path"].joinpath(gmx_xtc_name),
                 "md_log": ip["task_path"].joinpath(gmx_mdrun_log),
                 "cv_forces": ip["task_path"].joinpath(cv_force_out),
                 "mf_fig": ip["task_path"].joinpath(mf_fig),
-                
+                "mf_info": ip["task_path"].joinpath("mf_info.out")
+            })
+        return op
+    
+class MockedLabelStats(LabelStats):
+    @OP.exec_sign_check
+    def execute(
+            self,
+            ip : OPIO,
+    ) -> OPIO:
+        cv_forces_list = [_ for _ in ip["cv_forces"] if _ is not None]
+        task_path = Path("label_std")
+        task_path.mkdir(exist_ok=True, parents=True)
+        with set_directory(task_path):
+            with open(mf_std_fig, "w") as f:
+                f.write("mean force std fig")
+        op = OPIO({
+                "mf_std_fig":  task_path.joinpath(mf_std_fig),
+                "cv_forces": list(cv_forces_list)  
             })
         return op
     
@@ -323,6 +360,8 @@ def execute(
                 conf_tags[sel_gro_name.format(walker = int(ip['task_name']),idx=sel)] = f"{ip['task_name']}_{sel}"
                 save_txt(cv_init_label.format(walker = int(ip['task_name']),idx=sel), sel_data[ii])
                 cv_init_list.append(task_path.joinpath(cv_init_label.format(walker = int(ip['task_name']),idx=sel)))
+        with open("conf.json", "w") as f:
+            json.dump(conf_tags,f)
 
         op_out = OPIO(
             {
@@ -330,7 +369,7 @@ def execute(
                "selected_cv_init": cv_init_list,
                "model_devi": task_path.joinpath("cls_"+model_devi_name),
                "selected_indices": task_path.joinpath(sel_ndx_name),
-               "selected_conf_tags": conf_tags
+               "selected_conf_tags": task_path.joinpath("conf.json")
             }
         )
         return op_out
diff --git a/tests/op/test_prep_label.py b/tests/op/test_prep_label.py
index 8d770548..b344c1c0 100644
--- a/tests/op/test_prep_label.py
+++ b/tests/op/test_prep_label.py
@@ -14,14 +14,7 @@
 from pathlib import Path
 import shutil
 from rid.constants import (
-        explore_task_pattern, 
-        gmx_conf_name,
-        gmx_top_name,
-        gmx_mdp_name, 
-        plumed_input_name,
-        plumed_output_name,
-        sel_gro_name,
-        init_conf_name
+        plumed_input_name
     )
 class Test_CheckLabelInputs(unittest.TestCase):
     def setUp(self) -> None:
@@ -31,11 +24,12 @@ def test(self):
         op = CheckLabelInputs()
         data = Path(self.datapath)
         conf_path = data/"conf.gro"
+        conf_tags_path = data/"conf.json"
         
         op_in1 = OPIO(
             {
                 "confs": [conf_path],
-                "conf_tags": [{"conf.gro": "000"}]
+                "conf_tags": [conf_tags_path]
             }
         )
         op_in2 = OPIO(
diff --git a/tests/op/test_prep_rid.py b/tests/op/test_prep_rid.py
index 406b4f86..7e8015da 100644
--- a/tests/op/test_prep_rid.py
+++ b/tests/op/test_prep_rid.py
@@ -14,14 +14,7 @@
 from pathlib import Path
 import shutil
 from rid.constants import (
-        explore_task_pattern, 
-        gmx_conf_name,
-        gmx_top_name,
-        gmx_mdp_name, 
-        plumed_input_name,
-        plumed_output_name,
-        sel_gro_name,
-        init_conf_name
+        init_conf_gmx_name
     )
 
 class Test_PrepRid(unittest.TestCase):
@@ -34,7 +27,7 @@ def tearDown(self):
         data = Path(self.datapath)
         ii=data/self.confs
         for idx in range(self.numb_walkers):
-            os.remove(init_conf_name.format(idx=idx))
+            os.remove(init_conf_gmx_name.format(idx=idx))
     
     def test(self):
         op = PrepRiD()
diff --git a/tests/op/test_prep_select.py b/tests/op/test_prep_select.py
index 32dba677..e6e0b9cc 100644
--- a/tests/op/test_prep_select.py
+++ b/tests/op/test_prep_select.py
@@ -13,16 +13,6 @@
 from rid.utils import load_txt, save_txt, set_directory
 from pathlib import Path
 import shutil
-from rid.constants import (
-        explore_task_pattern, 
-        gmx_conf_name,
-        gmx_top_name,
-        gmx_mdp_name, 
-        plumed_input_name,
-        plumed_output_name,
-        sel_gro_name,
-        init_conf_name
-    )
 
 class Test_PrepSelect(unittest.TestCase):
     def setUp(self):
diff --git a/tests/op/test_run_label.py b/tests/op/test_run_label.py
index b1bd3c21..84ec3f4d 100644
--- a/tests/op/test_run_label.py
+++ b/tests/op/test_run_label.py
@@ -14,16 +14,6 @@
 from rid.utils import load_txt, save_txt, set_directory
 from pathlib import Path
 import shutil
-from rid.constants import (
-        explore_task_pattern, 
-        gmx_conf_name,
-        gmx_top_name,
-        gmx_mdp_name, 
-        plumed_input_name,
-        plumed_output_name,
-        sel_gro_name,
-        init_conf_name
-    )
 
 class Test_MockedRunLabel(unittest.TestCase):
     def setUp(self):
diff --git a/tests/op/test_run_select.py b/tests/op/test_run_select.py
index 975ee584..b6a2b4a9 100644
--- a/tests/op/test_run_select.py
+++ b/tests/op/test_run_select.py
@@ -14,16 +14,6 @@
 from rid.utils import load_txt, save_txt, set_directory
 from pathlib import Path
 import shutil
-from rid.constants import (
-        explore_task_pattern, 
-        gmx_conf_name,
-        gmx_top_name,
-        gmx_mdp_name, 
-        plumed_input_name,
-        plumed_output_name,
-        sel_gro_name,
-        init_conf_name
-    )
 
 class Test_MockedRunSelect(unittest.TestCase):
     def setUp(self):
@@ -53,9 +43,11 @@ def test(self, mocked_run):
         cluster_data = data/"cls_sel.out.npy"
         xtc_traj = data/"traj_comp.xtc"
         top_path = data/"topol.top"
-        model_path1 = ".."/data/"model_000.pb"
-        model_path2 = ".."/data/"model_001.pb"
-        model_path3 = ".."/data/"model_002.pb"
+        model_path1 = ".."/data/"models/model_000.pb"
+        model_path2 = ".."/data/"models/model_001.pb"
+        model_path3 = ".."/data/"models/model_002.pb"
+        gmx_config = {"type":"gmx","nsteps": 50,"method":"restrained", "output_freq": 1, "temperature": 300, "kappas": [500,500],
+                      "dt": 0.002, "output_mode": "both", "ntmpi": 1, "nt": 8, "max_warning": 0}
 
         op_in1 = OPIO(
             {
@@ -67,6 +59,7 @@ def test(self, mocked_run):
                 "trust_lvl_2": 0.03,
                 "xtc_traj": xtc_traj,
                 "topology": top_path,
+                "label_config": gmx_config,
                 "dt": 0.002,
                 "output_freq": 2500,
                 "slice_mode": "gmx"
@@ -82,6 +75,7 @@ def test(self, mocked_run):
                 "trust_lvl_2": 0.03,
                 "xtc_traj": xtc_traj,
                 "topology": top_path,
+                "label_config": gmx_config,
                 "dt": 0.002,
                 "output_freq": 2500,
                 "slice_mode": "mdtraj"
@@ -97,6 +91,7 @@ def test(self, mocked_run):
                 "trust_lvl_2": 0.03,
                 "xtc_traj": xtc_traj,
                 "topology": top_path,
+                "label_config": gmx_config,
                 "dt": 0.002,
                 "output_freq": 2500,
                 "slice_mode": "gmx"
@@ -112,6 +107,7 @@ def test(self, mocked_run):
                 "trust_lvl_2": 0.03,
                 "xtc_traj": xtc_traj,
                 "topology": top_path,
+                "label_config": gmx_config,
                 "dt": 0.002,
                 "output_freq": 2500,
                 "slice_mode": "others"
diff --git a/tests/test_label.py b/tests/test_label.py
index 43aef1c1..0a1b5dbf 100755
--- a/tests/test_label.py
+++ b/tests/test_label.py
@@ -39,9 +39,11 @@
     make_mocked_init_data,
     make_mocked_init_models,
     make_mocked_init_confs,
+    make_mocked_json,
     MockedCheckLabelInputs,
     MockedPrepLabel,
-    MockedRunLabel
+    MockedRunLabel,
+    MockedLabelStats
 )
 
 from dflow.python import upload_packages
@@ -140,7 +142,8 @@ def setUp(self):
         self.kappas = [500.]*3
         self.angular_mask = [0,0,0]
         self.tail = 0.9
-        self.conf_tags = [{"conf_000.gro":"000_000"},{"conf_001.gro":"001_001"},{"conf_002.gro":"002_002"}]
+        self.init_conf_tags = make_mocked_json("conf.json", {"conf_000.gro":"000_000", "conf_001.gro":"001_001", "conf_002.gro":"002_002"})
+        self.conf_tags = upload_artifact(self.init_conf_tags)
         self.label_config = label_config
         self.cv_config = cv_config
         self.data_out = "label"
@@ -149,6 +152,7 @@ def setUp(self):
 
     def tearDown(self):
         clear_files(self.init_models)
+        clear_files(self.init_conf_tags)
         clear_files(self.init_confs)
         clear_dir(self.center[0])
         clear_dir(self.data_out)
@@ -159,6 +163,7 @@ def test_label(self):
             MockedCheckLabelInputs,
             MockedPrepLabel,
             MockedRunLabel,
+            MockedLabelStats,
             prep_config = default_config,
             run_config = default_config
         )
@@ -169,10 +174,10 @@ def test_label(self):
                 "label_config" : self.label_config,
                 "cv_config" : self.cv_config,
                 "tail": self.tail,
-                "conf_tags": self.conf_tags,
                 "block_tag" : self.block_tag
             },
             artifacts = {
+                "conf_tags": self.conf_tags,
                 "models" : None,
                 "forcefield": None,
                 "topology" : None,
diff --git a/tests/test_select.py b/tests/test_select.py
index 8c778ab6..d5f2071b 100644
--- a/tests/test_select.py
+++ b/tests/test_select.py
@@ -74,6 +74,7 @@ def setUp(self):
         self.task_name = ["000","001","002"]
         self.block_tag = "000"
         self.data_out = "select"
+        self.label_config = {"type":"gmx"}
 
     def tearDown(self):
         clear_dir(self.plm_out[0])
@@ -93,6 +94,7 @@ def test_select(self):
             'select-step',
             template = steps,
             parameters = {
+                "label_config": self.label_config,
                 "trust_lvl_1" : self.trust_lvl_1,
                 "trust_lvl_2": self.trust_lvl_2,
                 "cluster_threshold": self.cluster_threshold,